{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import sympy\n",
    "sympy.init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define symbolic variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "t=sympy.Symbol('t')\n",
    "k=sympy.Symbol('k')\n",
    "x=sympy.Symbol('x')\n",
    "m=sympy.Symbol('m')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bessel function identities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAAAyCAYAAACNvqh/AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAPa0lEQVR4Ae2d/5XUthbHF84WQHgdLB0AqYBNByFUAHSQnPzHf3tCByQVJNABeRU8oIPHqyCc7WDf9+ORPOOxx5ZsaXx39uocj/xDupLuR9K1ZNlz7+bm5uy2uDdv3nxUXn+R/2VpniXjgWS81/Zc+9dL5d31+M7Gbg1wNqfPxvuzsozXajP3yxajnjQpCOP5Xv5iY0wuJQcj/Iu2z9rHOLubqQHpz9nM1F3taM6mtobnyy/JRrK8P5uPohOzJBcE57A57+TE6EFQ0Ff5vx/Koq5hVDEMT7V90/GjQ2HjeYX5oq0xyjo3GT7Gc3+rAekPnTubrUrM7DkbMyh6GZlio+ven/W0Vv/EFBdyUJON+RGyCv9KOriQj+E86HT9WtsPCvBV298HA+5dUJwPOvVBPobFXYYGpDNnk6GvYwZ1NsfUdl5aKWwUxvuzPLUuDp3ChURqsjFtkFXwC5X/N23PUUSie6xwPGtOdkoHY/9YPgbGXYIGnE2CklYK4mxWUnxCsjPYeH+WoNelQWZwIcnibEwbZBUYw3olZTHqnXQKdxkCJY+Qd4RilH+TDH+evKOUkV1nM6KclS85m5UBjCSfzMb7sxEtlr+UzIWka7Exa5BVYEarD+W/zdA9U9ZM9bDAIcspDlPXGH5G5O5GNOBsRpSz8iVnszKAkeRnsPH+bESfpS7N4ELSVdhYXtSFYTy4iOsADEbI7ehYiv5Rx99rY9TL61JThppR8keFY6ScNCpX+LvonI1d6s7mdNh4f3YclmbajMkRsozhz+KAEb1K5aE4hG/n9HWMMca90/ZT2Dg+6BQHY85rVQByN6ABZzOgFCOnnI0REAPZyGWj8N6fDeix9KlcLqRfk41Jg6wyv9bGyuepES36ia59fqx4zb58pqGp2Dkrr7kJ+DEoPcp2f6sBZ7PVhbU9Z2ONyDY/uWy8P9vqruZeLhfyUo3Nec2SzpEtQ8go90Jbzspqkmrm9OWjrE+S00xdy2fE+0RbqotT3jzDznl+nSr/1oZzNnbROZuTY+P9WWWkFtuMxRHyr+LAwqzcL3JhiBkJY3xfB2VrN88pHqNyjDJ3Tu66GnA2XX1YOnI2lmh08zKHjfdnXR3WOJrDhXxUY2PRIFPYOEpNgiAjyrQ0o2pekcKQ8pEPPonJaHuOIz4fI0Gmu60GnM1WF9b2nI01Itv8ZLHx/myruMp7WVzIS202pgyyCosBxbjyTliOY9EWLhryT5vDxkijRBaJ5bi/QmCAuZMGnI3dauBsTo6N92eVkVptM6YMshhEAxgNayoWnrfwXWqmm3ddfHXpX7snp/aDHGQh191GA87Gbk1wNqfFxvuz+jxNthlrBrlZyCCDGA1pKhZG1VcxcDCovFPMs2RGx+21GCbBZ5Q9d8o7QfytC+Js7CJzNqfFxvuz+jxNtpnz/XIHA8bHNOJ7vIxWWWTVWfWsY56zEoaRJMarxP8UPw2y5KU75QXldpzOLV0hzaKyeBfVkb3WgcrEzYWz2bwr7my2FdHbzVYXnb2V2wx5yWbj/dndtTVDBrkxZKoUN6pM/K1ez9iFGv9n8F8qDEZ5kZMM7grZckfHi9Idifxfrilfl9pyp9BHxM6/pHw4m436nE2oRqoT3m5GmtRabYYsGWPjbeYWtJmeQQ4VKU7V8mGNnlNF4x1d/nO4M2ruBcw7cRGCNxUnL2qV0HFhGLpYZJClJ0ZzrNrO/RRor2CS4Ww2MzLoxtls3i5AF95u0MKAW6nNkBNLfZr3Z9u6YYkLuWrZDBpkBcCA4DqrnVWxuRvn3a132i89ko1KKi2XcsxxMR+P5kTeixNHMXunZx06m+0sirPZdvqxvs6qVAUjxXxYYrNGm0Gllvo0i1zQ0RpsLHFBBy2bc44G3AvOyei2I0PtMxp5IZ/FUjXcwyA0Zq5GGskyVU6emxM+5is5buWAzsbZ7FaxWD+93exqpbu/RpshB2bYeH/WqRBmuJCrXTaHDPJjhWu/lKUITFE/kV/z61WMInHXG8/Mb8yXlQw5my0JZ7NZd4FGvN1s68X+3hpthjzE+mmJTczTvo7WOl6DTdSBJS5NfekZZBndy0CmGR3rmJW9/PtR7Tvw+K7wt5C+BQ9gDy1khDw4mw4JZ7NRh7ebTrXoHqzYZsiINTbeZjbVwxoXctWwub/JX+c3rqr+R5WZ/wVmZS+LkViUxGtOoy7EId7P7I8G7l580D00ccTNgaV8LWXzQExeafucqV1LOohZPxk24gEX2gzb+7Cl6jw1XNTbMXxLbGa3GXGgz6O98O9v9GewYUSX6qyxscQFHc5msw9AXHL6NGtcKE7DpjdC1oU4QuYZanyXF8PKtDULugZXXuv8mcITBkMeX8+51P5HbVHxBDvkmpGowl4fCjB0XuF5PSvZKfy95MCbgMkjZMl+pyhRf7vJxLINTfnzhbHU1epR9hw2dCQxfm6FjPl3NrtUu/tRt9lsJAZD3NaNUI/oYFIWRjmbLof9o0VcJIxXP5u3I+QzW/hvbd/tJ3LgOJuN0vD+LMHW7Oo7cKF/S3XZXBB8DDbnAyWgYH8r8fYVHe1TKZnCxsA+1vZlIB6nMNzP4jXiaMMgc6dZZcpbcnMNbMxecV95aTvVXeE6z8wCOog3OLuXc/Zns1HaMMP4T85y5GRoLKzScjZp7YZRGKOvuIiyuQHW8VhbG1P95LU7xGZ2m5ESr/YUyVQnI5lq7g5xQYdL2DQMpC9uuKrYliaBnZ9jsOlMWYfCkYXO604hT3H6mVFyzykuS8kZee0rh1EVip9yTUWXnNzR25TcpderNsDUzC1hk5rGSDhnM6KcAmy4kYvvIo6kNHjJ2QyqpRnRxNFxdn+GSHHlBvbLjnhuZmM/uHP64K5FNk2eDub4SBcKtBn4YCsao56ZbYtcKMK3jkHWiTi1HO/U23Kq8JzD2PI8ZchoDp0jPoV/yM6Euw7XU8JOiCp2mbzEfBUTOlPQEjYzk2yjRR04m1YlnZ1FbNSeftcWdYxgDDSzUl86qQwfxHjOpq+fRVyiOHHg+TGPEPj+QjtzGK+P+NbYnFp/xszSnFlHa1yoQg2bfYPcTGeqkIc6Ap6R4gZHyZtLvV8SetA72z/xTziVErYfu96Z/RF/vZTGJddgM57i9qqz2epiaK8YG7W9x0oAeU+GEho452wGlBJOFeESOn0exf2gfZ4jpzqLbE6iPxMH2B5czzQByCIXsvz1XAXDAP6hDZ9p5zOd448j/iO/vfvQfvNci+ty3DHScfBvShHwdXOl/4PcGKZ/dXvmUPxtiOPvkfdmeuP4STccSL8Em6XZdzZ7GlS9L85GMml/tDPe+U/VeWq4vRJUPUQ3q7SbGlzQFDy08VGkz/JZG3No0ELw6KyxWY0LCpHOSH9xfyY5tJOH8lPsCknvO2tcyF/D5lyFInOTq3wVjsrINuhQTpCFsvYr6/7xkIzYgBlRr+5UFhSEWw1e0OdiNptiLPp1NnvqK81G8mg3/GNaM80ajs/kT3U6zmaHjfRVqj+j/f9P2zPJjP1X7Aue6nw8t5N6b9cMG5WB8uBiGTZHR/wtxUZZpiyPJC8+z2/KFo4ZSE6NnM1wQf3Kb8tmaJU1YeY6ViWykKKprEqI/Q/ypzoV0oth6JgWuVDAzxLCFFOUmysz5iNOb+TG3w1PI2C7rS7qMOpkdjmcTV910gl65XEQBpmZJxzPkQ/eADchNj/OZkcZpXbFgRExHXfUL6Lpz2jHf3GQ4GLcRe1G+aDD9v4sKFz6wL40NoZTOka/vHJ7pf2UfrYIl5B2UTZFDbKU8VZb84EDZRZDxnuUL8n4lFM8VjQSLPXZ2ZhIGg6QaFBzXWxEEd5cOWcqV2+R3GxhMyMqD5SHTr7RjY65u2zfGR8Tq7DOZkxBy6/R2caG3UqT3uE16pzNqHqWXmR26lfpON6Uf6/j5McJBdl4f3aApHSMIY6L97A9vGY7OkIuyIVcFWVzfqCcs0+rsCl39YfkY/yYDlrqABSn0OfKiga5vRObK8hCPHFBt87GAoy9PIhN6ocm9mK2h95uWlWU26HjlrSl7b8EG+/PDmAVI1a956x8j5JKcEFWUTb3Y+6M+FT+aAiXZIm7ltG7pATh3A2fCTjg3G06JmdjsyZ4u7HJhVyVYOP9WXm+JbiQq6Jsio+QF+qNF/ib95xlCFOeBbTJKTyKYYqPqT8MB1856q0W1/lUh7zVp5pTM3uEcM7mCEqemYSzmam4I0Sbxcb7s+pkZnEhVzXZWDPInwIGjGHWCFdKwnjyOgLvp13Kj88Vgsh0T3Ex6Bh2oLnbaMDZ2K0JzubE2Hh/Vh2oyTZzv3qxMxJQJWQagZHxi4xo+0ExxMhZ4rghwGXdFGyinOavs7HL1dmcNBvvzyrgtdpmTBnkoHdeKYgGcQ4K4i6dao4P6v35cZeAs+nqw9KRs7FEo5uXJWy8P+vqsuTREi7kozgbiwaZ574PdAeTbZQVJz4/XjrVzLR3/Ewoine30YCzsVsTnM2JsfH+rDpQc23GnEFWJWR0y8h08gtVA7gaIx5kNJdDpR4IOnxK4XmvDTdnKf0m5on+Ohu7YJ3NSbLx/qwiVottxpxBDvpvvpsthTHizXGd5y2Kz0j3IkeAwvKuLl8X41m2u74GnE1fJ1bOOBsrJPr5mMPG+7O+HkufmcOFPFRhY9IgyxgyOsUg5vyrFEriO9jNc99gzC/kJy/wUljuSDHgSz6goein65yNXbbO5uTYeH9WGam1NmPSIAcGGEX+VSpnlEwcnj8z7fyT/LdBVqoXR8e+mGtcY85mXD9rXnU2a2p/PO1cNoT3/mxcpyWu5nIhzSps7t3c3JQoUBUZMqh845d3iyl8Vac0GB2zGOw77ft09YS2nc2Egla87GxWVP5E0sdio3S8P5tgsXv5WFxIc4yN5REyeX+p7ZUKwDRybceq6udKy41xmqadTZqe1gjlbNbQelqax2Lj/VkajxjqWFxI7yAb0wZZxpHnv1ehABSkilM6PNjnH438QyCJGnY2iYpaIZizWUHpiUkeg43S8P4skUcMdgwupDXFxrRBDgXgOTD/TUolK+4kl6kdPrU55zWr4vm5TQKlM2djFJizMQpG2arJRrK9P5uJviYXspTCxrxBDgXBWGI0eY2pmJM8/gweQ/+smNA7Jkg6dDZGmTsbo2CUrRpsJNP7s4XIa3AhS6lsboVBDgV6Iv91KBinFjnJYfX2H9qead+fGy/QpvTnbBbor2ZUZ1NTu8tkl2QjWd6fLcPRxi7JBaE5bP4Pa+jiR/hhDJUAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle K_{0}\\left(\\frac{k}{t}\\right) + K_{1}\\left(\\frac{k}{t}\\right) + K_{2}\\left(\\frac{k}{t}\\right) + K_{3}\\left(\\frac{k}{t}\\right) + K_{4}\\left(\\frac{k}{t}\\right)$"
      ],
      "text/plain": [
       "       ⎛   k⎞          ⎛   k⎞          ⎛   k⎞          ⎛   k⎞          ⎛   k⎞\n",
       "besselk⎜0, ─⎟ + besselk⎜1, ─⎟ + besselk⎜2, ─⎟ + besselk⎜3, ─⎟ + besselk⎜4, ─⎟\n",
       "       ⎝   t⎠          ⎝   t⎠          ⎝   t⎠          ⎝   t⎠          ⎝   t⎠"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test=sympy.Sum(sympy.besselk(m,k/t),(m,0,4)).doit()\n",
    "test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Substitution\n",
    "$$\n",
    "K_0(x) = K_2(x)-\\frac{2}{x} K_1(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAA1CAYAAACp3CmcAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAVFklEQVR4Ae2d67XdtBaFdzJOAeFQwYUOElJBQgeQVAB0AIN//Mu4dABUEEgHQAWEdABUQG46yJ2fo+X4bUmWt3X2WRrDW35IS0tz2lpLD3vfefv27cmDI+AIOAKOgCPgCJwPge++++5Llfa34t+GpYZr3+j8D9r/fni9e6zrpPmqe264f3d4wo8dAUfAEXAEHAFHYD8EZJg/k/QHikdGnlJ1/kdF19pecbwS/qv0vy6lueM9+iV4/Joj4Ag4Ao6AI1AOARnle5L2u+IHc1J17SNd+0vxnbk03fNKx+jAPcWTvf+rbuLa91UJvJZvFMd4OYvVkQzA/kXb59p/s5jYLxZFoBSPzmFRWk6leEEr58a5KYvARUnD7vywUqPHut709vUsYfQZAfhQ+wznj4LO/6gNx4B4ZM9uTI9eygPOr1RkVMvME5J1X1mRyxDKCJxMsZ5tAYHSPDqHC2AnXCrNC0U7NwkELCR1bpp7CWNnRu4T7b/mWNgkd/qU52vlfagN40nAoL7R+c+bo/ATcCcNtuGltqzygrgmkkzq8afiD7rnh/uh7D90/oU27BQ6/qMNW/W34lHQ+aZeinv1IOHdUeoKT4RKs2ihmJGnmpLHTcLN8yfHHvZFYA8encPtnO3BC1o5N87NdgRaI98sONM99ZU2hrxpuzGY9Hx7QecYwqZ3i1EdBZ3/XpsZQ+zKp53jbvrnOsDQ/iekSXYqusLCPvYmxo5Zve6r7BfacDYeKZ408kE2cj9TmnvhuI2qN/RSmrmHjxSbN9cq392hctro8f9P21/da0v7SguRAEnP3sNOCAjfVR6VxjncCf85sTG8kNe5mUNwv/POTYvtf7XXW1UubLAHGL+pdhsjiZF/rW0yKD+9ZALt/ygE7E+KS0/t0g4uLpxTmeiOsf5YG05Io6viRUdD18GDNE+09ULVhj5UGJLN++op3z2gkto+1Tk8nsmVjN303X3l46bBc4IED4UREK7cuKs8Kp1zWBj7JXGxvCDDuVlCsvw156aHKYabHvqwp0o7T+eA9qUbGjvAPds9Odi3HnPP6FKGNtqq3xRPOgEDOdGHkmfOxcuVTOhG+Tg3ODKNM6PjYT2nxCB7ZC+rNvRSGBKeqYJLwxXDygJmj7xhgpljjD2vKQxvppnkfjoBgVQencMEcDckTeWFopybDYAnZHVu3oOFQWeIfc5wD9tsDOWakX6KeMlsO4Xa597+VjFz8Sk2B1ExAb2W6mEycFTMhnVHJXqjGpZ4EDOazRqGXqjW0AtoetfXiidfF+jVIhworXlpLXlT6abOKS83BuTizXkohEAqj85hIeBXxKTygjjnZgXUQpedmz6QwoPhc4axhwHDfNK1V9oea/tFGwaSni8jtByzQG0qkLcdClc67A3z/3T49goPJXjOWemWiePSOCrSBx2pH/o97yaa2Sc9oxI95+dqJnENpzG4MYsWurriCTH8GwNmN5/tQzLz/PTs9/DorJzbFKfy6Bye5+5I5QWtnBvn5jwIrJSi9hlDjUFvDLOO6dwx3M0qeYw+9+pk0LVeh1DHOAM8D3u3+Rjfbg99Tr+e7tJvNBQ/mfHdSZMPNhj9JlRp6APwgPIs6BkbQWDbm5ccSMeLQhbDMYsOgK5zowAOpKeAq+QehggISx6gVB6dwyGQhY8zeUEL56YwF0Nxzs0Qkdlj5q1ZRD0c8cVItgZuJrcZ0n+Vn04dtoERgy8Vs2p9cdhf12nTWPBmbwDMFDM6fa0zezsTZuMoqw1VGnppx1yEvVLQKru0E8DHy2s+RABhIT3HvD7HFjNCgHPBkA/DHwaaTnnIQCCJR+cwA+G8LEm8UIRzkwd0Ri7nZgU03Yu06cx1T3XGcEYXDbWuk4bA6K85CnTuGB7/VttsfqXHxlh+DH5KIL31uCfzSX7Un88o3dwX8yblV2foA5AMO0yROAlOOGng0ytv9hXjLEAMXlTb018S0kkH6XYTrGTxy0MEMnl0DodAFj7O5AUtnJvCXAzFOTdDRMbHwoh2mbVb1itvE+kchhTbYQvZ2muDHWwCdqLt+GkfxwEbwbA/8/uTowLhPHPm1pEciF48pON4vZRCcucM+FK27jWT3zP4V90UlezjUeFpTQK9oCPEAyQN0kvlbwx7kDP7TeGhPKWnbPLiWbuhHwIUf5zDo3MYj29uyhxeKMu5yUU8Pp9zs4CV2mWM68eK206g9jHsJ8V05hpnVPttp077vZFZHZvDOuUM0KvnOjy0ZWj/JgWcHQK2sA132716dgC6JSpBLfJBNkaduRO8ttzA/A8f6WluolwhtzxfDo/O4f43TQ4vaOXcODf7IzBTQmjPHyoerorH+FvvFWe07SAqLdeGbThpCCMbo/Scw4ZMfl2OTBsDss0QR4uSXnwBcFiPufyTPfqqDL0qg3EGiClva65iJ+UjD0Dwzj09cQw14OQa+5+Vn0Dj5iERgYB7Eo/OYSLIGclzeKEY5yYD7MQszs08YMKGtp02nd45n8JtN52jU/cm5MbIYUxPOtfYBMWt4Q9pmiH3ifPhcvtnM/TqS4c/JPCTFKGhHkwlNPWKyAtWjEobJk2Wqgy9NDLDOvK2VirICkiC5Xv57vCdN6dKs/o7OgSQAMq8v+i8nrBBIIdH53D/myeHF7Rybpyb/RGYL8HejWd+frh1DRq9fZwB0jxR3Ey9KuZc9x37Uzju2QWdY+jeDPzXOuZVawxnqYDTgS44IbGBZ3borCzlZUTb7GCb7qrdq2OnmQcUELHei2ndDNkoX5d0rpmcDy1hQoyzkDsikFDMRSbN4dE53P9WyOEFrZwb52Z/BGZKULs+9bGcUepgN7hXeyHYhdU5d6XDURhODfRkbTmQfBYAYqPo1Y+McVe20uFgMDrNCAQ9dJyQ54rXjD6yf9DWC7X16FHSeuM9RVcO8JCeWZoAJoQxrIPX1l6zNBExgJb05iKKvJgkOTw6h/vTn8MLWjk3zs3+CNyOEpgWHjkjw6rLbvEWADYMx4Ap6dW/yFUanlM6pzb1rN13YdSjV2IM40Ntn4U0eB54FD2PSMfMmZAGRTDOq4oozWwISqKo9cJn0w4vKO8IOJ1rhm2GaROOm3/Akxxet1j0vhJkZiVV+Tgc5mnSWL/mWOcnvTudP4RDKqey4TCZR+W7aA7BhqB6HsJNLi9BZ+emYe/dj7D0tq+Dx23aFfe0xfS0GVJnwTY9bT68E2tv6G3/rs3ac+0uBgx3rP15orST35+ZMvSNwlKcF/fxKkYPeVDreYi/UBqM/dZgvefov5jdWuBKfhtZSAF6RWT6ZWELLiw+aXnQPjcXiw35C8PRTaBzR3FIBWvisQoOAcXCgdzUxAtwODd2U9T1zFTJzXuojt3T80tHNNZIj5RVft7BZwg/5gt8OBN0smPtK3q1dqJb+MjQc1GCbW76RTex7es6ix1eK+718u16ZmwNUXKPPrO8tWymR9T80JIw4WTeX/uBhqX0g2sYdTzINkgevXk4oGfxQXuhs6PrR3CIBjXxWAuHHWb8+QpgODfv74qanhm0KsbN+yr6XgeBL7RPr37SvnbStZ1MteeMkjKCMDeKiz2gQ2jcdcScTnd7R+8PMEyE3mtuFKYNw4NHsqZkIyDh5zqknVQ0QU6RpKqfeVGm1xa5kMSWE+Ai5b+YrYwjOKRsw+twHivi0Dix+AhuquEFEJwbuxWa+JK56VXUD9p7n3l3bOlSoJPJK3kE3gaYbFMlB0eRf/ibnT64QsJEeMo5ZWyHhbV/X6eeKs4etpgop3vKDKEZ2O61I/dNr6N0gAPeo5zDZU6/IzgEI9NnTt8jcDSdjih7qswjuDEMauIFbEyvKZyOOOfcvEe9Nm7ea3bD99SeM5fOp3wfa2vt7KBazOezoJzeOivu555dbPLi6Pqcoceot0MEoaAHintDyEpTMtgrcCw0qyUArHnbh+gkzOcIhKOTrrc8DRQ8gkNUqI3Hwzkc8MLhEdzUxgs4ODegUN8zg041coNeFxPUdi9O5Ya2fdXmKt1qmpGhVyYbVmy8DB2zSpghhslhg4Ko3ysoq5QonI7q9BInGAqGa/DkRuFADtGlNryq4vBAbmrjhXvFuQGF+p4ZdKqKGxTykI/A1By9rdrjlQH+q5dxfzwPFgLYK3ezJSoN8/hfauNvYVPCNYmV701KpjOkbfQ6QzkpRbAIj6GfuTmZbA4ls1mHoRju+ZpU85e9CcrVyGNNHG7hhmeQZ4sVu3y5C27uR3JTIy+ofhHcDDkQLyntn3MzBNCPiyIw6tFLuvXoWdZvhoQePfMELAiYXYQXGh3Lf5YehMqM+v9e6d0Epd/6N4AmahRLNnMqVv/udXuQp4ZYeN1ibni+K6PZD2XM/RezpTcdkjmUAAx8q2coj0Zr89sHptwwVhm3icNN3Ag7uG+G/BQz2sbq3ck3L4Y45xw7N81o5mrb18U28HK/e26P/cq4ob3nXkxp91lA9moJm9Q6Lsm6zdeuJirPDXrEf/UyVHQSsfQoo3v1Srub4Q7YNHqF/cVIurQGsptQ5xkJoTdmjlP3cvS+8tPgsIDDeoVzebM5lEB6jPQUbYFI4+TpmAWBiw9lUCaZR8m9NRwKoy3cPBsQzrx77P2ZzAtlOTdx/1NuvAgvHLnUac5L4IY2+4HhUCo+w/1XStWq5fSG7sNNisK91+pCDexVAHr1ewQz7td7CM+UiS6mV6aIMtnEDc7C6L+YdZ65+jYU4BBn5WUrMH3H8KqFx5o4tN581vMlbhn9edWhhHvCnsvO6cnd2nhByYvhhsqIm3uK7mszJ5nTMcG5iUHJ02Qj0DP0kmI9xdGNqpuYc3iqe/1X77+hFjwsNYVU77y47sKexmPtv5it3E0cqqwftVnDg0wMP8PFXQNjZU3FNfJ4OIcBqE3cGNjigvl5plP4QMbiyl3Lo7hGXlDvkrhhNCxn1M656dyovlsegaGhbxbb6Wada9SZgybs0avvGpd3pRz/i9PRDKsdpYq4oMfO4jumNNr/YWZf53jHcohbMQ4lGwcDeSlDckN9lP3QcDiHndoX4Ua8YEweaeMTyMzTx4TaeEHni+FGPMDt7PqlFYKcmxWA/PI2BK50g/Kw/aSNuBkG1jkMyx+KW+9U+7YgT5dO9CgwAhiaUh65GdRrCjg6qF7gQTj6IWSYF16Ynx+GxiELuhblUDIpE875fkIKBtXwWAOHe3DDTQAn2ni9kv88YE1Ncy9wbSZUwwv6Sd/Dn6+gw+bnRnJ4Vlg7k9sWOjczN62fLoMAhp5GfHXVt9LRqEy+t11GlXYIr3E2CsncIsb0sGG1LbKy8wr31dXupTmUPOrO9/SboeZwfFIc05BZGsMvu+4FMpoOh3EozIo8X5KDYfxH2yPtm1FHNuETbXauOTHxUxMvqHcx3KgucMP6GVsvwfEpHNNhWuvpOzcAdkOC+IRfOsM8d/zny2obfXTVplbdH6KTwGKhEWWnDBNP6hqI+FMXGdq0h2gy7cJJa4hy83dF0yCzVR+EF/VmWgBDz6gNgXn6KCdPeYrwKDk8TM4h6CsID3rw9Py69yOL+7ivfta2GJS3CC8UIlnOTQdtsNVh62jpmGeIETi+Z7763JNfm5JfZNtHvS4qBE6xLbRPWxYunw2XooZeFecGxyjQAPE6GR5uyn/10ojhJW0NTfkSQsOYG8zQtw9wriDhMFrcmCvrDPm4ea0hb4tTHeA1NpTg0Tkco83I27fi4t9w6aHilKmVErxQtHMTCBhG4gYDb4su+SbFr9rWevSIuUhuhvhc2DEdoWc3oU5XJZXUDc3NGtXzmykXo0ojsjXwoLFSfNWbXiiIRvQU6rSQ7LIuqb4lPr5SgkfncHBriRtw3eJ4luAFrZybATd2KI54CyL2TQjLRuzcdNGofF88m526EZ24u5XhycIzVpfTo9wSICHGi14qAxk3gsSlShx0rQSPzmF58krwglbOjXNTHoGbJRFnl+m0LZ3Js9X46mwlxRVk8x3JDUnwsBhexklg2J0vuY3eHohRQ/nIjxwaRg/pCGTx6BymA52YI4sXynBuEpFOT+7cpGN2ZA5sVNsR1PPB65WMAmM3WN9UlQNQlaEXOCxKAaCn2l5oiw7KB+i8ZgTgjxXjceUGSCQk6fAui/8K+ywencN9751cXtDKuXFu9kXg5kjXs4AxZ36++a6MjrE5BI5Z48SWM32jbPuE2obuqSUriM3Q5tQaA79lHpMybQ6SNQce8hDYwqNzmId5TK4tvCDfuYlBOS+Nc5OH27lzmX2iY9nsK6ZTiAOAzWh7+udWbK68Gg29fQXOwJzTfe48+bYCjYfWeGtzhfj5VQS28OgcrsKbnWALLxTq3GRDv5rRuVmFqIoEOLuMPPMsvA5G/qSYkUzegqmug1idoRdIGGmA4lWipKC8eFQfacueW5cM+wJdVUMvSUBUkDiXR+dwX/JyeUEr58a52ReBGyMdA4+N4psvfB3WvjdSbQWqM/QBqeZzu6FhSQEPAmiQ2h59hgxeD3yhfHhsHrYhkMOjc7gN85jcObwg17mJQXdbGudmG3675g72hM4kH0Ni8TejMHyGumpjX6WhF2j0pjG032pLCb35Q8lhCB5SooLS05CRfsu3AKLKug2JMnl0Dne+OTJ5QSvnxrnZGYHqxT8JGlpn0t6WaOyMnq3YP5k6a0WrNPQBAYwtf57DcHxsuFbCZn4k5OPrfCkL86w3X90cSywAFaZL5dE5PA+JqbyglXPj3JwHgXpLaZxd2ZXhiK/ZjA9rVL1aQy8g6dVjpFN69TRefHCHefYnir9XHBWUlt482xdRGTxRFAIZPDqHUchuS5TBCwU6N9tgj8rt3ETBdFQiOp7tZ2/FFQaf54K5enrz7TXtVxPuvH37thplhooIOOY9fte2+0pGlfWXyuFDB/7u/JCIjcfn4tE5TCPqXLyglXPj3KQh4KlLIlBtj55KqnGgR4+HtOurbiqHBTC8GuFGXkCUDufg0TlMZ+0cvKCVc+PcpCPgOUoiULWhp6JqJBh+55vCGOPiQXKbIXvFya/zFVfmggXuyaNzmH/j7MkLWjk3zk0+Ap6zFALVG3oqqsYCI/xYsX1qsEj9JY+pARyIR0UEupBFBPbg0TlchDzq4h68ULBzEwX/YiLnZhEevxiJwI0w9NRFN3zRjxNIHosqftL2SPvDFZQU6WEHBEry6ByWI6gkL2jl3Dg35RBwSVsR+D+D2qK0BLgCBgAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle K_{1}\\left(\\frac{k}{t}\\right) + 2 K_{2}\\left(\\frac{k}{t}\\right) + K_{3}\\left(\\frac{k}{t}\\right) + K_{4}\\left(\\frac{k}{t}\\right) - \\frac{2 t K_{1}\\left(\\frac{k}{t}\\right)}{k}$"
      ],
      "text/plain": [
       "                                                                             ⎛\n",
       "                                                                  2⋅t⋅besselk⎜\n",
       "       ⎛   k⎞            ⎛   k⎞          ⎛   k⎞          ⎛   k⎞              ⎝\n",
       "besselk⎜1, ─⎟ + 2⋅besselk⎜2, ─⎟ + besselk⎜3, ─⎟ + besselk⎜4, ─⎟ - ────────────\n",
       "       ⎝   t⎠            ⎝   t⎠          ⎝   t⎠          ⎝   t⎠           k   \n",
       "\n",
       "   k⎞\n",
       "1, ─⎟\n",
       "   t⎠\n",
       "─────\n",
       "     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAAoAAAAOCAYAAAAWo42rAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAz0lEQVQoFXWS4Q2CMBCFhQkIbqAjqBvgBsoGOkf/GUbQFXQENjAyAhtA2AC/V3tGG2hyeXdfH71LSzKO48KWc64KeYeuiQrWiiVmBLyoL+hDG2iGiO3J2zTAM5qZKbAB1UdX1d6IHolGIFpP6kKnm7EA9JFJpZ8PLdIwy4TnD+U6MQ9IM82tb+s5g/GlTpyazQzWrdOM1lL3Fi9jn3tktyZWsYvaTqzteu7A7YRxA2vU1RtJboAePZiZXG1L4iT2+9ba0E8xEPopdoTe3r/YGx/SQ0OZAIYmAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test2=test.subs(sympy.besselk(0,k/t),sympy.besselk(2,k/t)-\n",
    "             2*t/k*sympy.besselk(1,k/t))\n",
    "display(test2)\n",
    "sympy.N(test.subs(k,2).subs(t,3))-sympy.N(test2.subs(k,2).subs(t,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Substitution\n",
    "$$\n",
    "K_3(x) = K_1(x)+\\frac{4}{x} K_2(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAA1CAYAAACp3CmcAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAVSElEQVR4Ae2d65XetBaGJ1lTQDJUcKCDQCpI6ACSChI6gMU//mWRDoAKAukAqICEDsKpgJzpYM77eLQdX+SLZPmzZrK1lj9b0taW9L6yt272d+fq6urMnSPgCDgCjoAj4AicDoEffvjhuXL7R+c/hrmGuO8U/pOuXw7ju37FI/NNN2x4fXcY4H5HwBFwBBwBR8AR2A8BGeavpP1znUdGnlwV/rNOFzr+xr/gfpT873Myd3xEPwePxzkCjoAj4Ag4AuUQkFG+J21/6vz5lFbFfaq4dzrfmZLphkuO2YF7OkdH/+dd4dqvVQl6Ld/pvKaXM1sd6QDs33R8revLWWGPLIpAKR6dw6K0nJXihVI5N85NWQRulTbszk8LNXqs+Ga0r3sJo88MwCe6Zjp/5BT+sw46BpxH9uzGjOhVeMD5nYqMapkZIF0PlBS9TKGMwMlU68lmECjNo3M4A3ZCVGleyNq5SSBgRtS5GYMT2tY3Os+uTY9TNu3yW4U/1IHxxGFQL6Xr68YXfgLuyGAb3ujYPMiUToz2W53v6zzpQt5/SeC1DuwUZfyvDmzVPzqPnMKbeuncqweCN2KNPlSaTQvFjDyVlz5mBughvcXvbl8E9uDROdzO2R68UCrnxrnZjsCkBgZorGGPnNodU9iMbjGqI6fwlzrMGGJXvuz4u/Kv5MHQ/ifIbJ5Jli7szRo7xoge90B5v9ZBZ+ORzlEj30he6/1KMveCvz1VP6JXoVl7oOc2uZ5BbULlIP8LHe/l/4zwNU6yP0ruU52N/DXJXCYBAWG7yKNkaKDOYQKuW0XX8EIezs1WpNPTOzdxzIQLI1dG8n/revTMVhijcJ4j93WNgRw5hTNKZoCH0R9NhyuM5xV2BENfzEkfr7nRsYhuwiMjxdFBeafDOgTsql/VyZAcdULe0spb+Yg+VBgjPCKzKX3nR7JMvXypIHo8kyB2krSXSgfR9Jwg111hBIQrDXeRR8k5h4Wxn1O3lhd0ODdzSJaPc27imAoXDDTGO2rAQ6rGDtBm41qaUBsx93arKw2zATyr/tC5tJGn7Lg316fJ32Z9XvnTmaHDwnEmf3SGgriOQ/fIXt7tCNR4CQkvVMG56YphuQGzR95QYMKPsec1hXsT8R6cj0Aqj85hPtYpKVN5Qbdzk4JwvqxzE8fuqZ7RvdFqRAxDuWSkn5JOutpBoa5p29/rzFp8is1B1RpHuVgqmOuAoIeOitmw9wQEt2Y/AjMBzGr3XLWGXmAwur7QOfq6QK8WwSNZ66W15MXkYmFKS8OAXHpz7gohkMqjc1gI+AU1qbygzrlZALVQtHMTB1K4MGUf3a1O29Txmw4MJCNfZmjxkybmMOrtdLjkbGlxNI0fS5wZ9lDplow8qhlsNh0VlYsyskRB+V7pWHLIMyvRG7CeL6U6MB6Du9RzGxaPnhDTv2vAHKbFD8ns7Gdkv0ePLpbnbQ9L5dE5PE2LSOWFUjk3zs1pEBjkoucxxptne/S5rHAGd0y3sz6P0aetRp3iegNC+ekMcD9EdUeV5AVifLsj9KiWYdnlH03FRxNeB5p+8MLoN65KQx+AB5QXoZxrTxDYjualB9LpRaGL6ZjZDoDiaSiAA+kp4Erc3RABYckNlMqjczgEsrA/kxdK4dwU5mKozrkZItL62ZC9ZrSNgW8NXJu6f2GdgH+lk0EdtoHN2891Ztf65LS/4niefR/UYUxxzxQ+a1uuxZq3BPbuTFg5LkKezalKQ6+SsRZhrxR0yzt5HQhgOqaZ2pEfI4/Dz05EjjUzBHQumPJh+sNAU5C7DASSeHQOMxDOS5LEC1k4N3lAZ6Rybgagqe0xbR2dsh+I4qUzOmmogzwyOGYIbGmYwR35YMTn0tMxgKPG6drsy5q3vOgk2Ig7aOifpG/Vn89IbuqLeVH9d/vZHO9TBTDW9JReJJbGyGNU3lzrDGGASy+qHekv6DU5SHeXiUAmj85hJt5rk2XygnrnZi3ImXLOzRg4YYItYNC1OBKWDM965G0j21jhdQg2BjvRDvyCfp79rO0TP+UY9du9gAwdBF7Nnktjuhg4XpgndpaeO2uOWNoQZvp7Bv98JsFRUfSo6GktTb8My9esHyoQEt4ofWOwg57Zd/C7iiRP3qSl12a9va6IX69DIIdH53AdtlukcnghP+dmC+rr0jo3Y5ww3A/1TG5eMetEY1gxsISzk51p/cYAh+d3I6prOgkYWPObkY51BjDaxMPD1NItduGNjlodnR1cW2c8d/mpzAG0japTikY6en0YddZz1vSwpvTTeGhENDJ3eQjk8Ogc5mGdkiqHF/Q7Nyko58k6NwPc9Axm5M3/kfQOiWHILM7W7umMtgNEpWH5dvgMRwY3sjHkpXBsSPTrciSSzPBb8hh+OhptvshNOHSbIZ4QGQdLN5/MHdZjLHgdEh3RV2XoVRmMM0DEeltTFTtTOtIABO/cAzyGGnByjf2vSo/jxnOXiEDAPYlH5zAR5AzxHF7IxrnJADsxiXOTCNi1nRgaTYxcM8VvbVbnoQFu9m5Fwq0AP4ULRvWzLnCGvrUzxn9J9otZpYPIUA+WE5p6DaJjXuwgs9KX3ciqDL0KZoZ11NvqFjpy/SSEWTqbWqHSZ6r01LuUEVWNPCBxWO8vKueBkwjk8OgcTsJZLCKHFzJ3bopRMKnIuZmE5kOEnuV83tUGgo913X1XnpE9U/Xsr3qic7P0qjNh3Xfsz4K/ZxcUxtS9Gfhv5edV68aGfCjB9VUIRz7lD9HodFCWYQdlqL7rp10MOyvd+OE1nQ6zg23ceXtVx0WzDigg1vZerNTNlI3S9XoxijQ9n5hgwpnOQu6MQEI2t1I0h0fncP+mkMMLpXJunJv9EViRg57xzNhGXbAbtNWeC3Zhas29lZUcHQVbBmjDhxeSw/jzSl6TV/Cf6Wz2Zpik8SuepQZsFKP6kTHuJpIceVBXZgwYodOpeKXzktFHt81K6PLa3bWLSs4U0kbjKUWih/TCEgQwIYy1enptbZzJrDgDaLQ3tyLtxy6Sw6NzuH+ryeGFUjk3zs3+CNyAHGRPsAkYUmYWbIc+tub9yuKzLDzqjAzTSrdtMKRjwJI0HYtZI6947lMGp7b03KodjegljGHkIzNfBSl6HvQoej0i+VkHR4aCYJwXCyKZSRcKSUFne0UxBUo7Ak5hzbRNTH5lGN8MPpMepodme18r9WWLKX8al/U0eVjTqCbxlvwhHFJB5Q2HyTwq3a3mEGxwquch3OTyEsrs3DTsXf8IS3/2dfD4yC75HgvPN86tU5uYnGloha4v6CT8qcOe54PokRfDvdb+PJFs9PszMUNv6xq8uE+vYnSTh6K8Cue1XwUK4pMnGz03BnZS6nQRNrOQAnTx0gn/pgfZ5UHXTOOw2TD6d4cKP4pD6l8Tj1Vw2G0UB3JTEy9A4tx8aBjOzQcsqr7S/Xt/SwGVnu/WM4U/+wU+8pAM6/MMsi9X5knnIWqvR4YehVJsa9OvYxkons0O/Fdvb5Qfk00Is8aePKJPyCNF1Mqx5otHs3qFE4Txul77gYbZBP1IjHqvtyg9jObhgJFFtOEp/ggOKXlNPNbCIbi07iBuauIFLJybtkVUdc9Uy80HuG781TPVgFF91L52atcOMvXMYBYBGxKdvlc49oDlBLuvOmqm36PHMOFsd2PjITMdGB56JEuFbNIk/FwE2WhBE/QUEVX9rBdl5dqiF5I4chxcvAP7QWKmc+DDHuCD6PYNhlNySBkMr8N5rIjDGriphhfAcG56TcK56cFxuz2h7bPuji2dcwwyeSUPx9sA0Weq9GAD+M7A5HJ1dESvRE91nClhuzag6weE67x2bQEVKc4M2WVKohPIWrlOkFU0Czhg08cULlPlO4JDKmDlmSpvtJI7B1qZds5mtfojuDEMauIFwKxcq8HbWdC5+QBwbdx8KNkNv9LznLX0Cx2PdbR2dlAt1vPZUM5onR33U/cuNnl2dn3K0GPU2ymCkBHvC/amkCVT0tkrcGw0q8UBrPW2DymTMJ8iEI7OFN/yNCjgERxShNp4PJzDAS94j+CmNl7AwbkBhfruGcpUIzeU69Y4Pbtnl3LDs33R5kpuUWZk6JWIqWJc08uQn13CTDFEpw0QLOTuFdJTUg2djurKJU4wFEzX0JMbuQM5pCy14VUVhwdyUxsvtBXnBhTqu2coU1XcUCB3+QjcjSS1XXv2X73M+9PzYCOAvXIXSXYdJBn+xo+DLwstrUF09VzgUZrLbmAF1025KihHtwhswmPqZ2pNZiuHrP0/19F7haRbgJnrGnmsicNsbsRHs0dGZ+4vvvTV/J3yDBfdqBp5oXy3gpsu0FyLm5R7x7kZAuj+ogiMRvTSbiP6SzVWMyQYbNYJ2BAwuQlP8sjQQWjS6cz6A58RtIeboss66V71/72Wq+Sn/sfXRLLP0s2aiuHX1WM3cmyKhdctpqbnuzqa65AHrz3OpbEy5HD4oFOHk4wCVZePicNsbsQLBr5tQ6EtYFA2vxnSNK7Ij3PTzGYuPvu60AkzZkG5j3Z1HwM3qXXcFfAbrPzO1VX/GRuAZVd9zzjLz+5tHlKs1f8dq7PC/6fwR914XZPBZzrPTv0rHiPJDXVf11WM6lUO+2jOpgep9DATwoyIdZzkTXdKDz68Oz9n5M8UD+ZZHFqppIMyY1iS6i75qnhUeWrjMJsb1YW07bcT5Gf5hvpN3pOKa1xtvFAolenWcHONclMnnpF0kJltWTWocG4MPT/vhUBv6l4NzkYbvVeyQuY2Dc+ofuSUlocODXxo0DHaa3q3Ztyb0e8og2MCKIuV65gShFyFL4aXDlNr5ME84N6WTf5sDlsl2y4Mr1p4rInDrdwwmn+TSU9tvFCN28TNme69e6rTAx1Tu6ipc8w5NzFUPKwYAj1DL602ih81VDViwjDiU//VSyOPOTZ1rHno/xsST+mJ6T5F2LDjcoo8e3kIex4eD3Uebr7D+INv123hsKsn97pGHg/nMIC5iRvxv+W/sGvkBVhuBTeB3+fiKGfWzrkJAPppHwTOB2oxHGdqrNGpeUUxLcvInlH90OgoKOow8muMt/Vqo0oOCqTcQ0N60qKIC2ZK2HzHVDz4dx17IIYPlj047Oa5dF0bj4dz2AGsGDfinc4f+j7v6J+7rI0XynpruBEfcDG5f2mOGMU5NwsAefQ2BM7VQLnZftHBGaNypjAMy186t0ZE17YhDxF21POg4WV+65FPNVb0mgxpp5wZVDoGhzvVi3Ljpup1Hbv/L8so8ML6/NA1HbJQ1hIcDvXn+KvhsQYO9+BGOmkP3I8p/4VdDS80qtvETeCDj5+sec5R/aFzboaIuL8oAhh6DFm77julXXKM4CdH8TTyoIuHUGOAOrqG/k5Ue2k3SdPZaEOPu7By2LTaISURpoub4QLumzksVMGaeDycw9LcSB914r8OmmWA4D/T2XCfotHiDZMpuVOFWzkOu7+EWZFnnwBjUMD+GdvHhB9O8DNgWhrpOzcAdkOc+IRfBsNf6OA/Xxaf0UdXbTh1v7U8L6SADUeNYRcAXPO+tzXkSf2S4TUz4tdORc7pgoi3OtihvJj3hCJ7EOWm76q9lIfj1rtSPEqPczhoLcKENsnyDYaeGTUcG/QmO+CNhH4k7/eXgVH4DLZS2Q5m5IcnZuD4nvnife/cFCZkZ3WBU2wLNiZ3c+zOpeyrL2roVfGXOngli54sPXV6Os/6Wc76MKr0krY6OhjcbO83KDJD397AubqEx2hzY66uvdOprNQb49FgaFzq3C7jrChDCR6dwzHQPFisA9TGihv4WuNK8EI+zs0E2uKieQU2RPMs5DsiSyN6xJ2bCUwrDqazzeC2endeuoRq1Iuji5k8Mao8RLY6pjVtKSFX10MSqj7cgB+NC/XdwiFYleDRORy0OnET/UvigdictwQv6HduJlAWR3xFdPYb5hNJnZsJYGoMFs9mp27EIO5uZSCy8YzPfDJq2eIgYU0vei4PdNwIEucqcVBcCR6dw/LkleCFUjk3zk15BG6WRjq7fHl0cWmmhmqd11CIThlsvSP5QRJ6WExh0klg+pm/dh29PdDJa/JS6UiPHh6M7tIRyOLROUwHOjFFFi/k4dwkIp0u7tykY3ZkCmxUOxDU/cHrlcwCYzfYQ1NVB6AqQy9w2DAEQE91vNax2ikdoPOuOYA/1pkeV66DRFxSGa6T+K+wz+LROdy37eTyQqmcG+dmXwRujnbdCxhz1ueb75rIj83B4WcfDUfO8o2S7eNqm7qnlr/qMEObU2sMPOtdW5ytQX5U6/NbAIuk3cKjcxgBtFDQFl4ognNTiIiIGucmAkqFQWafGFg21zozKKQDgM1oR/q1lL1GQ890O+v0BmYqVqTbCjQ9tOFX6FLL8bHLb+HROdyv9WzhhVI5N87NfgjcDM10dpl55l54H4z8mc7MZPIRq+oGiNUZeoGEkQaoxQ/ASKbnlJYe1ac6stfWpcO+QFfV1EuvojfAk8ujc7gvubm8UCrnxrnZF4Ebox0Dj43imy98Hda+aVFtBaoz9AEp3sN/Hh4sKeBBAA+kdkSfoYNXy/jIDz02d9sQyOHROdyG+ZrUObyg17lZg+42GedmG367pg72hMEkH0Ni8zczZG91XbWxr9LQCzRG0xja73WkuN76ofQwBQ8pq5zkeZAhv/U98lX53XahTB6dw50bRiYvlMq5cW52RqB69U9CCW0waW9LNHZG99a3NdagSkMfgMLY8uc5TMevdRcSbNZHQjr+rz1lY56N5qtbY1kLQIVyqTw6h6chMZUXSuXcODenQaDeXJrOruzKcMbXbMYnNRa9WkMvIBnVY6RTRvU8vNjIxzr7E51f6rzKSZbRPMezVQlcaBUCGTw6h6uQ3SaUwQsZOjfbYF+V2rlZBdNRQgw828/eiisMPvcFa/WM5ts4XVfj7lxdXVVTmGFBBBzrHn/q2H0no/J6p3z40IG/Oz8kYqP/VDw6h2lEnYoXSuXcODdpCLh0SQSqHdFTST0cGNHTQ9r1VTflwwYYXo1wIy8gSrtT8OgcprN2Cl4olXPj3KQj4ClKIlC1oaeiekgw/c43hTHGxZ30NlP2Oie/zle8MLdY4Z48Oof5DWdPXiiVc+Pc5CPgKUshUL2hp6J6WGCEH+tsnxosUn/pY2mADsSjIgpdySwCe/DoHM5CvipyD17I2LlZBf+skHMzC49HrkTgRhh66qIGX/TjBNLHpopfdDzS9XAHJVm62wGBkjw6h+UIKskLpXJunJtyCLimrQj8H/fGqF+LmVp2AAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle K_{0}\\left(\\frac{k}{t}\\right) + 2 K_{1}\\left(\\frac{k}{t}\\right) + K_{2}\\left(\\frac{k}{t}\\right) + K_{4}\\left(\\frac{k}{t}\\right) + \\frac{4 t K_{2}\\left(\\frac{k}{t}\\right)}{k}$"
      ],
      "text/plain": [
       "                                                                             ⎛\n",
       "                                                                  4⋅t⋅besselk⎜\n",
       "       ⎛   k⎞            ⎛   k⎞          ⎛   k⎞          ⎛   k⎞              ⎝\n",
       "besselk⎜0, ─⎟ + 2⋅besselk⎜1, ─⎟ + besselk⎜2, ─⎟ + besselk⎜4, ─⎟ + ────────────\n",
       "       ⎝   t⎠            ⎝   t⎠          ⎝   t⎠          ⎝   t⎠           k   \n",
       "\n",
       "   k⎞\n",
       "2, ─⎟\n",
       "   t⎠\n",
       "─────\n",
       "     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAAoAAAAOCAYAAAAWo42rAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAz0lEQVQoFXWS4Q2CMBCFhQkIbqAjqBvgBsoGOkf/GUbQFXQENjAyAhtA2AC/V3tGG2hyeXdfH71LSzKO48KWc64KeYeuiQrWiiVmBLyoL+hDG2iGiO3J2zTAM5qZKbAB1UdX1d6IHolGIFpP6kKnm7EA9JFJpZ8PLdIwy4TnD+U6MQ9IM82tb+s5g/GlTpyazQzWrdOM1lL3Fi9jn3tktyZWsYvaTqzteu7A7YRxA2vU1RtJboAePZiZXG1L4iT2+9ba0E8xEPopdoTe3r/YGx/SQ0OZAIYmAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test3=test.subs(sympy.besselk(3,k/t),sympy.besselk(1,k/t)+\n",
    "             4*t/k*sympy.besselk(2,k/t))\n",
    "display(test3)\n",
    "sympy.N(test.subs(k,2).subs(t,3))-sympy.N(test3.subs(k,2).subs(t,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Non-degenerate expansion:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pressure for particles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI4AAAA0CAYAAABcrAAbAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJrElEQVR4Ae2c7XEUORCGx5QDAC4DkwEfERxkAEcEmAyg+AX/KMgAiICDDOAiwJABXASAM+DeR6uWpR2NRrN762E9UpWs75am+1V3SzPrg1+/fnVzhWfPnr3S3EdK7wytQW231fZD6ZehPruo13yPRPe10tNd0N93mpdmfoAPmv/b0BoktOtqA1gbg0ZjXwzRp17t0P+geDnup/JLlYtj4/5Ly88NHDTN5wLTX0iArwvtNU0JIOIBog0wHTjj+ij/Tn2Oo3LLeg4czsyJm5ofcBwpvctalGenk1LX0zSqR9CMu6b4XRFTckP1D5VOChoD/S9Ks5pF9R8V36nPtuCdtK596Dy3xgEEBFKEE/s61H1VXA83JUz6vlW84/O79EMGNdb6wn73snh1rPhVEf9tqzCqcTQJjHviZ0ELEB6ofithabyBht1u9GLgXFX9DyaLg8bZ7sdpxkfqVPfY+igPvVjYt1WHE24B8/PRChUp69iboGfDtH7LPaPqcPbhT0+Txw+oPq8Uixp8FDgiiCkJRCCqOvwSTMU2AcFjlhw90b2h/FWl5iyfqEyfJKj9kSLjAJkDjMrsJAcopQFEDFR5lAn0K4QeeAt9Z23Ss2LuMdu2uZL1qJ6Nf1np2MZB5hwY4o2c0KoxVQglFiCI5SRiGiMhOKEA8BAqQHmv+JdioKl6dsUtxfVwS23sKrTNfZ8HZJODxl5XRG2jVZ4o7/wsI0S78sXdaX3nTrVWZxmUhk2eWRNydKBRP2TIJkSeSVAdMkEzD5q0w2REvsBCNhJMntyqVosKD6h8oiWicaAeUxN2iPL3ovatsqIFKIjOIc8Qe6g+YZ2Z9t+pCic+Nsm5taFB4Cmah02BZvpXZdvAKq6C6jBr+EPZu6xR4DDQiPkURmJDd74T/aLZFcxnJmxtOaPFjXwxzYf2GRNEmFz92Z1oSNNagP1U9QnQVUbA9GFdbMjHqtuKlxoPEDg0DJoWtRPQOJ8U0bRo+U7pn4pDvOX53ygmz8C4g/jmWARQd/gvnFZ6xFQHSnlw7OhGAtHYyYF1ned8LHDTOTWOq3iAnvUBVQ9o7ivagUDZ7YJoImBAOqS5O7UBLk6ppgjQMkXAqh08/FS8onwi73WNAyKZ4IdiEvzE2MOtQSNak99zaEyynvMorM+p8kFpXrWzsQhuN6+yZ3/Vjm/G65PeDj7rtVEOujXahnspzC9yRgFcUx5fp6ckWIXqASPgwv80wNHUrQOHydktCbpUBkyoVLc4X+6UZieEcCloXFEApbG/eRsCIbhrglV2pb2Uf6LY8yWsz6apeGlgHfNDkZ2tK1YMuB6Dmkpt0AXoReDw4Mlu0cIADaoQ4NgixyZT9/MPfq2scy6HFhPUaf7YmYdn91VXEs42zEJmvc2eIYjZcbLVWrgtJ6Kp3mb6xlWYNzROEg41mIlhNIQBCY4TauyTUk4b+Dy0kYagtqxwVE9fdtd33xlbjzBPfXmXCc+Su23e5ZwxbUAS/AY9M4LZ6HVITHQkf0vto7zVWhJTpnKtueR58DETPxPgsDuwfThtt5WuT3BF9VVBYxEc2ume8o6BSmEmQNrVjhNpt8s51QBm9yyaN9Gc9Nll8M/OFE7bqMx68Ak3MucQqgxs1Nj0VA6r7ma0USpOpow8jIYDmNAQ1VdlxSgIY0NxvmI6gMlsaxWtTTppzpeKc9672Ib7rnVw84qWRdseK72rOAhktZmW5tHhI6H21HVVfXcJTtNmzBNCDBwEPPhwYcRwBk1DQKWx2wh/KHLhFGy+q93BH898e8gdzDBKEv4ROInYhSIaB3OFxi3xdpvXOqMaR+upOsWqX+7QYhpHj3AWHHA8001jnLVOyzngRUybNnr73sx/AhmtgV2enAK2Jz9KAZOMmQzzKo/TyqbBBcB3jDVxTJD1xi9fHeBGxth4NkuiDazBUtHJAcKax1KjnQDokh/ldosmCJpBeZA8NXArOVdwQtG62eHhOc5jMZrT8U9z5UwyICCgdYYCvpkD/VCHGesNBwA0hEOfS/wbMQJHGbs5tEMCgShDf0xTL0BPsaSqe2OmVog+82dPelNpbdDf/JseYLUutBBrgwfJycTmUX3QUr6O50Bb1fAf2iZcI1lMRZcTMgcYxo6Fosah0RERMRbBbWLNouNJ2Vm281w9tBTxfabSiunuQ56N1hV4Zv5fSeu45xSN68pAj89MagJa/mZNR/qIPvLFbNaAhiG4MPhtWY3DUZlbTdR8p9ScO4pVQWN4iwpQAJDdpVCeSwtUrXvTTjyrxr5RdBsNOqqL77+oos75K67QdbywBRic/nqCUx1Cov+U1zpsSvhMTISr+lxgc0/ZyAC4p0mTl5y5WVrd+XDAgybcevtyp7QHsPUVqQ8vIjE9PQFbX7UBSjYx2gyA0fet6osgUjtmDaWSmFPzcdTWwlwc8ELFnAEcNBIBIddemv6tvvhZJeAAQOijcZ4rHfU51Qdtynqgn4QGnIQdsxXY1QiJNAQJrtbMA7p/FGuABhAGARYmX2V4R/Ve6+iZwGaq1ji1r0UJF/8K0zOoSdSGtuG+qOo1kvrhq2a/zbJ7nH3lV1v3GQceKDt2agvaRqDAmTazeEbF59TGQWnwM5AGnB7L9rNCgsac4LtwKhsKvDuzS1pAlnW8RQNHGmd78HTdgDPE4j2sl6AxU3xgjknKBXwhvvpDm2DWer6LH4SvVPzsovk4nlMtmcaBpnGm8av19hw4ePr0adUr98axxoGYA81Uxdxo+WoONFNVzarWMeZAA07MjZav5kB75VDBKh1beR1gl2tHfkjtN8EVM+xflwacOplt801w3Qx71mvRpgpNosjH9D8V7RuinAiP1R5fqnE7+3/8q5fcXHtRt3TgnAoQfI7A1XvpjTFvqU/2QqLntMhmqlaM5mXf8yGeC1zJR0zqB5BqvwkeIrvX9YvWOEguMkEljROErP6AbMo3wWHsRcosHjgSJqYKkzX0wi/IW32OVJj6TXAYf5EyzVStfpkRtI3AgTbhh/wcwfnU0gHKg4YyQOt8mTT7aQJ9LnJYNHAkdMCB6eFzA0AAaAiU+YyTyK830DTUbfpNsIZerLB0U2VHbH405/JK+aYFQMUnLQBEu4GJlCO600bKLy4s+iWnBI8W4YNsPlw6Ubn4U5HFoaPwwE3jrDQLPzrjR3KD3+AWeLjIpsUCRyDBHOG78J0u9zL8SuBzA4+4UBEWCxzxBhNFsBOV3QwDpk4Asv/xQ7GFNQ4sGTgcq/kHiusOrh2vs/95Y41/iy0uGTiYqvCawQMIJxlfB20T2haLjsKD/wfwx52Saf/9qgAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\frac{t^{2} e^{\\frac{k \\left(x + 1\\right)}{t}} K_{2}\\left(\\frac{k}{t}\\right)}{k^{2}}$"
      ],
      "text/plain": [
       "    k⋅(x + 1)              \n",
       "    ─────────              \n",
       " 2      t            ⎛   k⎞\n",
       "t ⋅ℯ         ⋅besselk⎜2, ─⎟\n",
       "                     ⎝   t⎠\n",
       "───────────────────────────\n",
       "              2            \n",
       "             k             "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Pk=t**2/k/k*sympy.exp(k*(x+1)/t)*sympy.besselk(2,k/t)\n",
    "display(Pk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pressure for particles and antiparticles together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANwAAAA8CAYAAAANHtQDAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAPEElEQVR4Ae2d7ZEctRaGx1sbwLJEcE0GYEdgk8EFR2CTAZR/2f9ckAE4AgwZABHs4gzgRoBxBr7vI+uo1GqpW5rp3undkap69XV0dHQ+JU3P7L0PHz7selqOAy9fvvxR2H5T/utyWLeLSeu8EHV/6Hmk8vvtUroNys62QcbdoMIb2+WpGBtS80b2VMU/Vcb4eprgQDe4Cea0dEnZngn+gfKvWsbdBVit+a3W8b2e13dhPWuuoRvcAtyVwn0uNCjcowXQORTC+V8995fCtw8ezf9t7TjB/iTYC+iuHXOKcN3glpE6nv0nKdvkGUb9z/T8pWdSkYETvvfK/16GvDIWzYGRQNefKZTaftDDmbQ2Ed1fa0zfWhY41g2uwJjaZikXHp0I92pujGCJApd62IJlk1fWL5X/ngVoaBSOSWNRP3Rj3BhIyUh+ERzRezYJDofzRk/fWha4dV5ob24WsxHYcz/QtkJPvRCa8d2iASgjkQBlm0yCgS9ElCljAt+koUxOMuw0OQxbfU10YPhvlRe3gdCq53s90D27RuGDfqL4fT2rR2i/lIMz0Yrj+Vv5SDa+7zv1/6jyD6XJ1Ef/N6V+2hczOOFCKGEyJlcb25TP9NzJpDUSIVDq2ejmGfBYuROoxjIORf9UZYRp6bHqgY/WqDbbhsJPlJqxrqy+tRUbmqF99qMOaNEDPGuI16XqNpPohZdfKGcHMkq064HnxZ2JH4QN8JHQlyMkvuGs1LFHO+cAhGIJAvFyKOVdTUR0IkCN54cHCAKBYGzwBQHDN+o7n49wqR1j+1U53vVSj3naByrfBH+vNE9RidSXpl/UQMTYfBJPL0Tkc+UjJ2fEqw/5zO1MdoLD8bEFN+doKEK+pMFB8HXAfBoFPCPKVZvMIX0uoWBAGBcfGFuEQrDvMsjeRzAoyG8ehi37bNTJ4Gttgk7nFCoHEuFQUPiz9YT85rbwg52J1vWtHgLKKKkdJ/qNcuQ0SufWIgAYalsAPCeC/07tc2HUofATubL/gwGyvagaHw9coyw6iAQc5lnnG9WLHq1mfo8PULdFnBsjeOZFCGwDP1Pd8UZ5zB/6Ue5BEky81UH4jvZ4rMooTWoUfC5oxmk48cAxPmufytGF6iT8rA0nQlScdAiCIxo81GPGCT9xMIPPM1XHMICBPzj2at0UbDYJJ/yCR3PRO7cz+Z/GsdMwZxnPgSzQtcEaAHAG5ydmcJhYZSyYtwdGN2ZqQzE4n9E3mlBtKDfM+ULPJpJoQrHZp/+lnG3ZoQmntBO+0foLiJ2XFDzejzIKhOHFlwvggrfZJFj4GqJdPFblkQNR2+R5IjtJvhHFrF2nYcBwHI+sIZeLRncJoZx3DDHUoIMJ/M++vuRFHAGmxvkgrys9bmcCHaIz3pnQFCdwZi+azjwUxjUQmBBCDN4ExUgTBCCEkefTONrBh3IzfmupVXFK9KMYcXQqwVm785K+EvMt8F38Ah/8Cwl+6vnXNzxRHtNvO5IAv1IBmnBULQl4HMRs0voMLhsN1e/Og8q/0rOkToE33QEM6NV8rB0nyM6EAONoVV6UvaeR/q/1DNKZr2FAXOWCOE62Fx8ogQBQHrzRYPGqA0eohzA8McqSjo3x3+YyPIuVf24t8NYplHiCMLiOR+DmuW08fA08U5k52AKz9UI5uKThDEH9IIMTDuSDc+Tyx5U9XlUHCXlnjWEANaw43gifGdOwd1iDl6SB8mss50DoY82t8zuEpT8RXdclGN8ObcyPYyT4uACkepCRh0sz8Oa3lOrAsAiXAwOKMKSGCBEDBngC2LticMZkiDxIKSIatlaEJ9UGJ56gtCGpPhKG70TBBnwTbIiC6kNWiyThhf5J+QjGyd7DtsxrvEExi9HAI3xCrjnC2lRGh54on6TPj98ncw5T+Es6bzindiZTtBHhRxHuHKyatCT8ED4FA4EIHgHARAwUa79Szj6cMx19g1eE1Bcri7rHSTCMQ9HixYM3GHUEY1sbQjwXAKmQoBM84CQhtNFZUm2sgXMm6aEeW4drqPxjtFSCz4OxHj1snYg4prTzA/MQh44HqzmA/AzlVpv7sgwSetCzYJRaN5GfI8ms7gQM7QVkHutbCQN65PRQ9LArKe1M0vGshwjNE+ZxBpdCUhcQTEApnRWrjmKjDCjpY+Wpx/5E7c1JeJgDI+UwbAvDaDAmtlcsEFowbraqJsidylwKYPi86QFjOKimdOWEBixrsAM7+Lkg4qo+4FdbNgkGmknvPmbL/hV+LlZYyys9QVits4CndUwMr/EovpND3F5Thm49gMLrYhIMsiY5x6k6W2WMfFYODDogQdes/ERPqk+l4JSSYrjRFYzPpTMrZHIU3D5sjbshICCIO/YsMw+GFaKZ6iYkUzajJRUCzgDFxGC4EeOK18aq6hLb3DQBF26nVLb1mPBT+LSezpH2H1wXTVPblYPxVyLg7Bh2EJVjUrBP04akbgr9j+ZCljhBZEOEtx1IMuRjVf3ujOfH4Rh4amVzKSymX1n8BzYabuYJ6TyUooKIRkkxgpw1o5SxcUQj24qeORiLizQ2Wu3gt2hHP17iyvotFxzRj6rb66uMV/lXOUrCATznMNRc/IC+RVjgMaZSLibRs/fX6v36irjX7ijNr/Z7lXPP8dScHBHR9IAIR3TlMmdK1zDQEMVVRm/ZLXHcmEvQhb4Uk/BVyU1wOV5kcY8MToNZ6KVy8zyBILVBJMo/uE0KAO0FcJH++Zhl/xrMlHJjlCTOaggJz4jQEAjvwQWhqG2JZMycUyY3l+bPCWQJOu4CDmTHUSXeceDscZps+zkyvC0slNfiiGoWhZ2hzowxVOjTpVVyufAcIjfDbbripjiLJ9IEKCofxobIpjKh3ZTeeSPVbYE7lauULp4nKtsWccojGczUPAjI0aicW1LWALMwNIRiBhlN3Ys3xIHixZLkYtEt58AxHhIOtJSQ73Wp88jtpq8YdkjB4LxSPlSenh0wQrPSwflNsPSZMQaktQWNhxi814PcGPDroR+4XMQ1gXHGw6iIziFpLF6TLUkWfwBsLxgTzYu1Y1hhhNbLxc/e8liBJFAar3LoTabBgRuQ1kEbzhYduLD2OFd7+qVfDBDni87MJXBn8ZYGCm8Lf003zHYcWmdwQoSQUFoOobziFR61cWNmTAOJizhqg1iiX83iBFpMRFPwpMaChzPcj1T+WjAYVZyA4YbSBPZc5ZSJ1K2fscYIyvsmY2I61774Dh7n1832y3YEB+M8BIHoMMM3XuXQ4bB3gjU5pzCcyUhTUc4BCAe6Ab7RR0AOYPyHO4FqRyz8yLqFv6yfc+n7eGo7wxHSARgovQeMmUH0wxgdnHI75HrQ9kw48Ej/0UjOWzDYtiDMY8bN5QiMZLtoC4Be6mZMtBPR2EIqc4kbMvAyBwJBcOQwAgfzSg948IwkDJYoH7bUH5uHf9VvjJzaCg8HrV8j2seyWn/G6RlQUJLJy1XEO9pf6yGH9zu1IYsr5UGfVHbnMfqVeLMGueH8Rw5FbeABvuV1QnhFgOEZ0Kj2XGrlL/pquhnw3eu/Sxl40VSQkLgNe6fctkVN45cC1vzmMPDuKA5C/lntRzU+zQ89GNInKtcotEDbk18/jtc5TV/fKR8ZZopdMLyjyksGI8MwWPXtxV+NQz8IGuEyCJwW4Qx/z+s5wGF99OpO/fBlICVQFAuFwwPzQfnUNfoyk9ZheSgwdhZrGxu7ItZPBCRheOk9hOvI/HmjNhzmlME181e0XAgn9IB/kM4GtV5p4QAezG1JWgatCIuAi4qz4rwl1NCzdpRFBjgacns4UtQaOcaaO0apeZRa+Isjti8YDxD1CDdgR1MF74XAEPhRo4oUDBrsXKniJhI02dl4FYK07r1eJzRiNJ67Afe6ovKiDNXXyl8ibPao0SOccb8xlxDwokSULGMb0R0KHryv6CLqUj9a8grK/FuKuCV+PFUHl2lTqZq/WjsRM1z4pUi7waUcaatzKXD0c5xo4Lb0ypOO8sxeGHjYtTKcENHj2HTMrk804jg5+3LLWUpV/BUOLli4hAm3rSnCbnApR9rqbCuJKGw5jpnY2vJ2Dd6VG8raM8xaNEPHq7WQL41X/GI7yRewS3Ks5S9bycmPlPrHAgdKT0LCM/LOX+0HrgfOuO3h4gNfr+HzzIPOV9te5f7U9Qi3P+/cSCkWXo03ZY56bjpwGUsOZ0vLuainDAe6wWWYskcTRsfbEyedfHTjs7fijd9JM0iL7wa3gAZIwdzbBMprP9NZYNZtodDa74siotvkGWZbVN88Nd3gluM5L1jzxgOKd1JJa77Qgnkfl5/J2PzN5DGF0w1uIe5L0bgZ5Dq85Wv+C81+dDR8PMJL4n0rOSOKbnAzDGrp9t6dLdXJRDmtmeiGsQ1e0m3h2ynB9o8FTknafa1H50CPcEcXQSfglDhw78WLF1W/THRKTOlr7RxYiwN9S7kWZzvezoEMB/qWMsOU3tQ5sBYHusGtxdmOt3Mgw4HzTFtvukUc8NfyvOFBso8j+AD62N8Y+EhR/zvgQDe4ATtuZYXPwMI3q1XmqyT83MCWflHsVjJ2DaL7lnINri6Ak8ilh/8OxP9KsJ8OzGF+pv74e1x8Xah/eyHHqQ20dYPbgBByJMiI+I0SXhXj3cSpnyogul3ncPS27XGgbym3J5OUIr5nV/z2tIwyfaUKA+QrMmv/YlZKZ69XcKBHuAomHQtERmNbxakIF8gTPMbZ8nPfYWwv3AwHusHdDJ/3nYUtJVvL2RtHwdwXbOvPfe9LVx+3Jwf6lnJPxt3QMCJciG4yKqIXv2jMG/p8984Zojc26hjoztfJ+3fTYMiGUje4DQkjJkXGglGxReSaH+PB2Eh27c/VP/+uichG274/962hPd0UB/qW8qY43T5POL/JqFxZOV/wxBDjm0sMj377qW/yZ4J10U/lnjbEgf7y8oaEEZMigyFq8SOz/EDRter91jFm0C0t9wi3XcERtYhk/N4l/xet/wzfdmVVTVk3uGpW3RygjIttI2czfoKbz9X4zRD+3W03OjHiNqducNuUnv2/AruhtDdJMMKdDI9fN+7pFnKgG9w2hcb1Pv8MI734sGt+/pVyT7eQA93gtik0tpThdS5veFyecJYjuoW+bZLfqSpx4P9P7os5F+V//wAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\frac{2 t^{2} \\cosh{\\left(\\frac{k \\left(x + 1\\right)}{t} \\right)} K_{2}\\left(\\frac{k}{t}\\right)}{k^{2}}$"
      ],
      "text/plain": [
       "   2     ⎛k⋅(x + 1)⎞        ⎛   k⎞\n",
       "2⋅t ⋅cosh⎜─────────⎟⋅besselk⎜2, ─⎟\n",
       "         ⎝    t    ⎠        ⎝   t⎠\n",
       "──────────────────────────────────\n",
       "                 2                \n",
       "                k                 "
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pka=(2*t**2/k/k*sympy.cosh(k*(x+1)/t)*\n",
    "     sympy.besselk(2,k/t))\n",
    "Pka"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show that the pressure with antiparticles is correct:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAAoAAAAOCAYAAAAWo42rAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAz0lEQVQoFXWS4Q2CMBCFhQkIbqAjqBvgBsoGOkf/GUbQFXQENjAyAhtA2AC/V3tGG2hyeXdfH71LSzKO48KWc64KeYeuiQrWiiVmBLyoL+hDG2iGiO3J2zTAM5qZKbAB1UdX1d6IHolGIFpP6kKnm7EA9JFJpZ8PLdIwy4TnD+U6MQ9IM82tb+s5g/GlTpyazQzWrdOM1lL3Fi9jn3tktyZWsYvaTqzteu7A7YRxA2vU1RtJboAePZiZXG1L4iT2+9ba0E8xEPopdoTe3r/YGx/SQ0OZAIYmAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pka2=(Pk+t**2/k/k*sympy.exp(k*(-x-1)/t)*\n",
    "      sympy.besselk(2,k/t))\n",
    "Pka3=Pka2.subs(sympy.exp(k*(x+1)/t),2*sympy.cosh(k*(x+1)/t)-\n",
    "               sympy.exp(-k*(x+1)/t))\n",
    "sympy.simplify(Pka3-Pka)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number density for particles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIQAAAA0CAYAAABLjpDSAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJpUlEQVR4Ae2c6ZUVRRSA33gmAB0ygAxYIhAzQIkAyACOv+AfRzJQI0DNQI2AgQzACBjJAL+vpm9RvffrmZ6l7XtOTW23bi137ep+c/D58+fdEvDixYufoXuT/Ls++vTdp++E/F0fzhLtzPcUur+Qf1qC/nWm+dWCi/8T2h/66MOM2/QpMLOFgbE/9dG3nX7p/0n6usSj/or64NgS//9UXlIgtAxvBw7zJxjzy0D/lK4ao8sB0FbgktCV7UX5d3AeF/WtyAkcLngKd6Et02+SP3AeymqmuW0ty0C7DHTcLdJHkib9Du1PyPcCxkj/HXmnJaD9L9Lv4JxVKPda11VHXtJCyFzB3EMvYwnb3pOacBcmifua9F1VXtLP91qY5sKucp1zekx6TzI2OhMsYiFYWAiD2vmIukwtBeKI+klz5eCFthpsGoPsaHsWeJSlVzLxPm0GrwG6gb+iMiF3HdcC2Jfu7UPX/mgzQPZsWla3uTnwfib1WtxFBIJFyFDdg8x6ywLukB+RR5B5TF2cGtD/lOQ4hScJAnWlPwkKeRYOB1If3Jw4I9ASyhH8S+lmn7pcXWcoTG0dtOuCvyafogy6cQPtUkEzvaVchjGAzFIA/iD9QAqrsaNdSb5HasI9+tQErcPDqqzw7A2MvU3ShGoFfqSc4pggZD/lUY0K/MvKWacW0fX3ajX9KlcSBvB8slKxtBgtoF2eaEk73cvBUvcQrZU0GlhQrwlsoPZWoTHbQpxlbO+CFuhgnSqHDOy0Dk5Jn8HxG5LKp6ArHP+QtCphlal+AdqN4ezXnWdYymXkCfoKLES/pyTrFzsX3Te2aK9tpmgfLDKf1qKMPcbw1SYtWlgZD/wTdL4vB1KXMeK4Li3bM9pmWyHG6goMtDvNO30BWggFQquoUOzIvyUNnav7/5VU20PNQkBA8+TdgRH+EDFQzgeckzSLsXNXMHdOxnmtqwDrEltAu8LwkBSBdAtnnwboyTQFrxY7lTToU2jU9rAgWs1RIQRHXv9L+oZyPv+mhVDSnOCENAsgvvddOGNmzXWWQc05qR8M0aM/YqCkgU1c+nWBXsPXNK6Jt2ddmlOsg3cqT0jyTyt1i7KxRK9S06egKTjGdyFMu6ZAOLkakCWG+l7A2MGD3YvY1UL2sIX0OHxaTKZZTfuRpGb2MiDwp+bQCgHU9QyBPIs1lYpsENprWSqC0laAs0A0nzLcdKcGVAQuLeOAlPjJfn+BheoKdqwhP9pRlmk+ARgrnJswVGuXF1OUU4FMPGMNary3s1qW16Qx0NXcLZEOGezESpOEdRcGJpqdN+TpqpnyjnJoglfKgn7Ug/iUasv/cZ1dt5vLz3w6g8z3wBOwbw991rV6RWIsuwfC6NmyjppLob6Py3I/xnA5jlMglHh9kAHRffLaBLTtaJMZauf3lNOhkCftoG3MLIFyNmAuo3yFNq2T+oVaMeZz/0KyDtV6fM4/b6uQJqn+qIAnZcMC5aCvIUh8PSwmURBSY9G2Y/Mi66MMWsp+Dyl8Vznk3MvM+4rk/EOXM+c+b0EwlOQja/CmT8uohXxM/oDUK6D0hWWVnGcpTHkKOQJvSYFzHWGBnCtBKRAyuGtj4bc1K3G7dQNcrz+zPz0lt8zf6lBj8ctMMkzVsxGMzMONaiF0GwaUXedGcwIFKAsyZc/TR/vOR9fTIenvqIWA1qQnOvD6Av2wEHnaJBAMcPKwBLmzKiRBAScOotl/EXXXcOxErEOtzFHxRUzOHLpH3VWel7IBnwqhmzXuKq1nuSzXW750S4I0MsbxKkDW3JJglKHRx+hAGcuDfhaMr6oRSQOYIGs8ZYUkwFuwy4R02KxJjcxrvIgFMWc6G+bqco8yV9BK9IHWIQlzH8IltgePFb4Eh1Veix84BANM/ZeMMNdFtEA80pC5bI2Z08AcriGb3Tk0zjAm4oeWILIurYZr8xxypF7ORXu2KlW7+9C69FmUGC7dYFi0jebQ1R0Z/Dt+DHothB2JAIRchM/8sWC1ILQkTSAOSV8YOKl9pX9Ujl1xHs1teg7CkJVICNC4TUF6fg4wBlrlu2NIZT/05Z3ua4owONQwwbioZSF8dPSmTZO8I8/xAmVfQikACkbcA1i/LI11iYuC+2WCX0lJOZyMtq67mRQP2A/4ok6G+zTUYghtHr74rTeMtHWByuY5mzLDuhCLNhV3HyVVMGuWr/ZyqyC8Fc/xBCph8FE1KVFV35G3BKecln5fPmn+a0wrcSzTr7BJW+uj8Ij/mvZB4aBf96IhyG4tYgjaN1jiBCpm6VYUCC2IIPOmXOj9Bp4xzJhAKFjS10K8JB+N68DR+rke58iwCUQ+isUKaqGHb54BhkxxuQrS36QpwiNtGTwoPCJV4FvOP1hHzR1tLiOO54rmMMzYRfM/qPX0ax287/hmylbAMx5sffcS9xBTaGw4l3MCj5h29AkGnGwdYLaBaLin1qrp8+Gh83X9JhCt47paDTBPk25c4BPKEHgVHheIClBnwAodA1AD1fwkWRLdBKI8jStahnm6C3+Io1voA+MNv5RS+3UxtdigGGQ80vuKfIshipPairvdZiE2KaidwMHz588nvUKtjdoqqz2BzWWslrXzNra5jHnnttpRm0CslrXzNrYJxLxzW+2oTSBWy9p5G1vtyy0uZnyh5HsAPzLxJ3ZjH7XOO8GVjVqthfCmjuSrY69wp74BXBl799/OagWiOApf8nR9IFugbMU4gVULBBYi7v43CxEcH8lXLRDsXZeh6+h70TNyPP+/7tUGlRUrtRDZOiAYfnPoj2gNOP3kTGERR7diu98eHJH8mukjSYHyF2r+lFCchyTBsr/IyrRT6wr+rNZCwCyZLuNS/EBdYRB8TeznYybBr4b8NkB8v7T2NxYKi22+JpbxvlI+qtp9fSwNn2BWB6sVCDiV4wcYmcrkflcg49OTR9UeAedN+2grPxzRWgi2l9bAdumsDg5Xt6MvG0rxA1WF4TgYSu6n6emHMpS9n9BtyFwFovnhqxZGyJ+pn1bT+MFP3Cu8a5et3UJoCWS+P54J5mYmKQxVJSxIaQXsUqj8jyyBV6End9PEjb5rna9SIAqN91tEtV5/73/UbQlFxT2DxS6NV1BqjIeGbVoU44gddcurgVUKBNyJgDGYeVxxTLcgE59W9ci6GC+jxY8YI3AVMK2GP9i1P+aK/mudr1Ug+kx9fIl8I7hWMVXmv462Kk+Mpj+EKrrFjTZdUTO+CLxrma9VIGTay+AITPtE2cdFGah1yH2U1fKun+dLo4vZfg7vU4d0ktsgXw38B0sHhvvD1KIfAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\frac{t e^{\\frac{k \\left(x + 1\\right)}{t}} K_{2}\\left(\\frac{k}{t}\\right)}{k m}$"
      ],
      "text/plain": [
       "   k⋅(x + 1)              \n",
       "   ─────────              \n",
       "       t            ⎛   k⎞\n",
       "t⋅ℯ         ⋅besselk⎜2, ─⎟\n",
       "                    ⎝   t⎠\n",
       "──────────────────────────\n",
       "           k⋅m            "
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nk=sympy.simplify(sympy.diff(Pk,x))/m\n",
    "nk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show that $n_k/P_k = k/t/m$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABoAAAArCAYAAABmbJjGAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACi0lEQVRYCe2X3XFTMRBGrzMpIEAHoQN+OiAdMKECSAfk0X7LJB0AFTChg9ABSUpIBwR3YM6RJY2si23pwXlgvDPrXa1W+61WmuvVZLFYDC00m82O8LuG38CPjF+2rEs+B0nZJgk8h0/we4B/bvOv55uBioWv0G+KcZPaBcSO3sWoO9+RpbOE86ZtFE6Hhd6iuqO8GwDfM34Le1HONyXQXDqCGCyfTwTBNHyBTyM7/ic1A7E6nw8gQUf+wG4CW29iT+nC+UTAW0BCCZH32F7DG6l3R2Zu0DMALGMzNQER1PIcwxfoZ0i/EHc9YE1ABPWwpXTjbpfDAD4A+DmO14pWIM/nnoDzKpKllF4sxfrfViBLd5HCRMBzxp6Vu8lzyaeWk9avd72wd9y6o964I/890KgkrYZ96VorNfJ7stJNptNpW781yrHPsP8y9NWr8H6yy/D/AfV0QUXFN6v8GfpHeQefoId/4V2Vzr7vGH6EA+0KKDxv4l/+ToHckV1spvBlANkJG0Ib9o/wc9j+7Tdsdjf4XMH6fIAl9UtsqWM1hms8H3Xtdk2/XJuAXOBrwAdWnkQfsLnwD2yQB8YpsC+Jb4yfITMx1n6NnGQjygEG0dMLzgM8wnZVOLk7SXtqIB1rN4maQg9YGw8x2LD7uHKRQGZekiWSvi5F/rUHt8GvycRXzkcHd2SpJB0GxmXWmkKGhZ82yTZ5xReflGyqUHD0p7zeHvK6DOuAJmVQH2Empy6Nkk1zJZBOdcB1GVpee3Evh+VOj4CV82HOi+H8ckfR2aDfNRYUAjC/kgDz+iab/Xc6Py9I+ORgC0kiQ5XSjkQ1u7p0OqcgqJku0byFNvihfHHGxl/7J+QpMt/ev1uE19AECqmMAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\frac{k}{m t}$"
      ],
      "text/plain": [
       " k \n",
       "───\n",
       "m⋅t"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.simplify(nk/Pk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute entropy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUEAAAAzCAYAAAAJk7pPAAAACXBIWXMAAA7EAAAOxAGVKw4bAAANZklEQVR4Ae2d75XcNBfGJzlbAIQO8nYAoQKSDiCpIEsHycmn5Nse6CCkghzoAKiATTp46SDLdhCen9bXyLLkf2N7PDP3nqORLEtX0mPp8ZVkj+99/vx55+IIOAKOwL4IvHnz5q10PJT/pKRL5x7r3I38j6U0c8errBfS+Yv825zu+7lIj3MEHAFHYAICvyvP36V8IqGvdQ6SnESAyvdTSTfxOo/u3+W+iNPp+GcdF/M6CcZoedgRcAT2QQAL8EOHgp9ESL90nO871SC3OLH0QrCBZOP4KPyr0lxGx3Xwog55wBFwBByB/RB4pOwQ3UP536NKYawwfOIaFqDiIC3y/E/ukxzT1W8U/6P8UaI86P4oP2vxKf4PuV+VpkXCbgmOgtoTOwKOQAcCkBqCD9nEa4PE/V8ulkciJtK9l3tShbPrdnGmPcJZS9ItwT0Q9ayOgCNwh4AIzAgQS+y5jiGzmAQf6PjmLvXdb0V6HLBZwnriTnEv8RGF0RUT12PFsfliwhT3DzsY4FOHljgJtiDxCEfAEZiAAETG1BeS+iBy+kb+A/m2UXKtY9LUonMv5MgDWQby0/GlXJiyyq8JkUw6fis3eqpM3koaJGyRxemwCmMO7+IIOAJHjsBKY5l1PUgK0vtN7qmcWYc7xbNm961cLN8qns0KrMBnVRiyHC3K+7Ucj8Jg7b1SOKxJmiLOK0wd2sJzgql7/fr1C7nv03g/bmPlmDgmW+8DGsuPGdOHrqfqcEldptZDed8ukfceSmOpGBSGbpiiSRrY+wc5M29h/liwIjE9X0pPnn3j1FFY6VkDYBeHXSMequQOc3BRPc6uzVNAP3ec1H4sjh/lZ6dtS+LTdb1ULlbSrfzW7mhXvrnPVfX4Tb5NkwcXoTzsPBd5qaRIebAK/5af5aIGCSoRBPSnfObzvaJ0/yjRtfx4ATTkUxxERuFseWcLDwkLP8rD80boznamQrbFo1Wfs2vzFFDPFSe1mx1QHtXASCjKUvhIL2OYscNua4NodEz8d/JvixVb4QR1XLMOfeWlGyPsxsS7L0VIpBhrD8DDrk4mIXogQUhsCpFxR72S24ycY5ungH9MOKmuzGZ402BvC0k6sLZ6ZWF8QntUCWZiqTAmGeNTxmOqqz5We5rTyfpMOaA85ZMLnMmVp7h7FJWS4FOdGAoQYCOlLWrbWEmfDbrL1fGrOvTp7si96Km+ep1im6cAekw4cSPH7SXqs9y0sbBwfbIYPiqYWRlTv1Y9FMf7s//IsUzVOt9X6dJ56QpkUjq/9fj7VkE1JMyb7XiAH6bAylea6kKmAJ2uFw5QHS4k6xezXaghhQ5Ic45tHgBLK8k54sTu5lBrckl8INiuMccUmZ3bTYuwxDofNCvdtyEXkQIuTMmqi5LVQcBupVfFv1D8Ozl81h8a6xKKGyIN3dIBQbO9js5Z72JDKhOladTL4tduc1WevQcJLs/lsEKfVXX6S2m6BkKVbDFvEzgt1rpEsbBmGjxmwM6Kj8pHH0YH44N+wOMirMnTD8JrawqbMGYZ64MIW/nR+UqO19oQNirXGIO0afQskgqOlYsowyOFB11IAYPpDzg7heN1EAYkF4HnhToXhsmbE+VDL/pDXXQMASIcs7CLG3QBlW42UT221GZ2ycKyhXzw56ZDh6dzghPrPgchQZW/JZwEw7Ki9tLfmbUMutkvgY90Qmy8G8tY4a2KYGkWWg6xQDC9Ij2koz/9oHCY8cnn+r6SG71LqzyDRGXQp+nfoU06XrQvxyQI+eQWU3MVNxB5FKBx8XXMhXgnfyfXICsdU8ZTOfKVdqBNNwCEsHy21AGfslrWp+LWEKvX4DarzrSXDoMwWBB7peju6O7XdPe2WTqxACE5k1sFwBxrEHkgF58PkSv+WFsOitOK7aWdYwhhSXwgv9LylEHCGLe+aHEtX20iDZuetC/WSf1Lm6EtPVMiVN7PcpQbbvRTdIzJcxElZvAwoIYIYGfvfqo4hIVFiDXIIy4BQPmQmHUAyKEkQbdOkpb8gfQqPSXiDLqUBr1/ynXpD2mjn/ouF8XlgqPbLCW1xYZC1Y+7KpZs+uzjmDaDSXzjARMeyQjXTn6nBa7zS2KkqtTruXEdid+p7Gzf0KklcApl2o/KBnvrfxaNT7+nbrkB1/moi/JwQ0LvGFmyH9G+PquJ6zJkfFi7eJzFZntfKS//17eoISL91C/0Z/mLy8XEEvrAtrk86YwE8elUNr0tFU0eLhSDm2cMIdP4TlTKt1M6gOskymLm/hOj2yyVl6pT/JI3FhpxrNnEbRrc5iQftcayviIwRJR/SYyowiZwSrFQu3Mkt1M8/ZFF+HTtLFXROFZ6LCUIokX2jYTtg0XwUT0gDrPe2qX+FzPU2An1HIvLf8XsFaLsazSofMZLY0a5l+ZM5vtR3I3CANkpqhAWHdJlEnO3Q+KBfhfT8SvddiGvFKbTsrjLy9hWZkfu5U5F5Y9tM20IF7NUu33arLx0FjCr7/7oq3SWilwsXuXaddoUTgs1GMLhzSpucrVTHBg8ruIayxJL4kOZcjuVUVtpCufGM3GM9SHy15BEC6QJvKH6Y2nX7VmgnKDyIlLMHY0L2ycBbCXKVk4VpxNwd8Xqy6bpKACrBrF8RiDUC332rxMh0Yo/k9qs+qZ3MEgxfX1ncJuljw7MjYHpIxihL12W4OXxMWtUUjGbbAKn2VrToajC3/ppnVLxvFHE2m5uWWJJfDA8aqND5TMGGdN1nMLIA7kh1itpmP62BN1y9Y23lWDPCOmmbPr2KnIRlQJYrOX1CY9hpAMv5FHlucgMUnR9FyLH/YQLKT23STa7aNmLkqRd4nDvNqtNdnNIp+tj2gy+OHaBW3f0Cv9D3b3BfSs4UZdDCdcFl5Ml8anJreobTPFzREX/S4kxV1es2AYRVXqJb1i4uczHFBeT4HtVHALLigDgHBeXwbyrjgmaYK3dyLH7mQPf0nX56L+yBNID2b7UMTtFrDPW5yzNkr7KnKXN0gM2dBzWOG+TOo9pM5YH1iVEuJMucGGh2jZcbhSeij0qJ4nK3BpOk9qxT6bqGthMKkyHpY/HltjpXAMfxgnr5+H5UcottIe+0yC3XDrl5+0Sllbot7bGz3Fv3py+Lcelf6BAY4fulk5ql0DETGc6l+6QTtK39UxqJwMDyy10nup4J9+s2603YZX6VbisjpPKnbQxsgooUSFz4FPp4KZ5FmMvgq8zGFuCJIT1Gawnx/Y0bm2pOh1WGoM7WNAKgy13bZcKgQPjhGWeWuebujYz4kO/Y4y7RAg0LEHiBTg7e0yzZrVUqgsJAWCOQwiY658UXzLbdfq4RW1jkZzpbkMUf9QvnDcaM8OB49QN4hz4SAczEqbL9uRGd6FndDZHggxatvwdrDPqCN7U00ZA43kR4+YUUGuRII0SYNw12AY/WSvtFC6et8ERGIKAxjFvfPC2zqyzuyFlH0OaLAkeQ8W9jo6AI+AIzIHA/TmUuA5HwBFwBI4VASfBY71yXm9HwBGYBYF7+ozd6O8DzFKyK3EEHAFHYAMI+JrgBi6CV8ERcAQOh4BPhw+HvZfsCDgCG0DASXADF8Gr4Ag4AodDIH1t7nA18ZJPGgE9o8ZD+K+qRj6s/NynBk4aB2/c9hBwEtzeNTnVGg39C/1Tbb+3a6MI+HR4oxfmGKqFdSfHv5LwQW/7u6VS1S+VhvfGTXiRn/+8sz+WsHj3HYFVEXASXBXu0ypMBMb/PfKOOa9jtf5lOWktf55xncT5oSNwcAR8OnzwS3ASFcCau+pqichyyKcGulT4OUdgEQTcElwE1vNRKnKzKW6fJViDojyQJn9mmn5qoE7jAUdgLQScBNdC+nTLYTrMtHjQH5Mq3UOlZz0w96mB00XJW7ZZBHw6vNlLczQVwxKsrUCRHBYeH+zikRj+Ubsmx4oAiQv/VVkd7+SHv3iSjy6sRPI/l3sgx1riJznysAnDNztI80wOIWxf3wsR/uMIjEHASXAMWp62gYDICKKDhPiEwE7HECBiH376oHBYC9Q5LMC+Tw08UTpIkj8AfSfHh4rCpwjko4ddaEiVz5ZaPGXyIaMv5VwcgdEI+HR4NGSeIUIAyw3hO7shLJ+v3UGO6Y4xhEgafHOXFantqvyQHwJh8vhN/Ke+WIUI8bXlqWPiKc/FEZiEgFuCk2DzTBUCYT1QYcjt2shJPt+1bWx6KK7PUiM/a4sQGiTINDgWLE4k3WUe+h3du9z+6wgkCLglmADih6MQgPyw+CAiPs5lRDVKCYkhwCoTOjmOrT2iINyPUTrikKdyadpwwn8cgSEIOAkOQcnTtBAQGZnFdqUwVhvrch8UnkyEVSFseGBJpgI5NshOZRFHPWxNkrCLIzAKASfBUXB54ggBLDDEiMneBmEquxNB8XGfKZIjOyNcWzM0vZAv1iEbJZRrdbLz7jsCvQg4CfZC5AkKCJSmp+FxF+X5qpCvGF0RGYT3PkkUyE3njXDtNGktjul4ul5o6dx3BIoIOAkWofETPQhAQFeWRgR0qzCPrUBGWIH1OUszwMeaw6pLp8OUlSM4Hrpmt5jywpRYvosjMAqBfwEtMM1/vC4bfgAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\frac{\\left(B_{1} k - B_{2} k x - B_{2} k + 4 B_{2} t\\right) e^{\\frac{k \\left(x + 1\\right)}{t}}}{k^{2} m}$"
      ],
      "text/plain": [
       "                                 k⋅(x + 1)\n",
       "                                 ─────────\n",
       "                                     t    \n",
       "(B₁⋅k - B₂⋅k⋅x - B₂⋅k + 4⋅B₂⋅t)⋅ℯ         \n",
       "──────────────────────────────────────────\n",
       "                    2                     \n",
       "                   k ⋅m                   "
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sk=sympy.simplify(sympy.diff(Pk,t))/m\n",
    "sk2=sk.subs(sympy.besselk(0,k/t),sympy.besselk(2,k/t)-\n",
    "             2*t/k*sympy.besselk(1,k/t))\n",
    "sk3=sympy.simplify(sk2.subs(sympy.besselk(2,k/t),sympy.Symbol('B2')).subs(sympy.besselk(1,k/t),sympy.Symbol('B1')))\n",
    "sk3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute terms for energy density:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAA0CAYAAAC6jF4zAAAACXBIWXMAAA7EAAAOxAGVKw4bAAARb0lEQVR4Ae2d7ZEdtRKGj10bgOFmYDLwRwQ2GQCOADsDKP+Cf65LBkAEBmcANwIMGQAReHEGvu8jTw+aWUnT8+Xdc053lXZmpFZL/bakljQ6s7fevXt3CAoEAoFAIBA4LgS+/fbb71Xju7p+Wqu50h4r7VLXP2o8W8errK8k8wdd39Zk364lRHwgEAgEAoHAjUbgF9Xur1oNNfDfUxqOaZHTUb7/1mQTr3Rk/6JwJ+fT83d6buYNx5MjFveBQCAQCBwPAqx0fm9U979yAj800qeSBg4lZ5ZcnFpybHl8dv+zeJ5mz4Pbi8FTPAQCgUAgEAgcCwIPVFGcy11dP6PSume1wZW4wUpHcTgK8nyi8EaBrbD7in+m6yxSHmT/oWtxZaP4XxV+Fk/R8cWKZxbcfmaB/lThTwX2O4NOBIGw64kYMlPjiG2KI4G4MsDn73qI+1MhpwfSFb6XCp9299X3MHnGhffVFdNmKx6Mp8r9peuvcyupPFTweZcPTw19qfg9QXlfyoq/LZ2Vxss1ZgODWUepOPF9rzB71lGSFXHrEZAtqm1ZaWdrV+l+lP2UFrGFTSXD1U+9fGtaqsowp8MYY2Nl7ng+VvxlXoby2OqDAwe8Hzoo7mvj0T2ycmfxWHEcYDBi+2zO+E4dirSJ41FlWOaxZDPFioU1Ilku9gOv7lGWvUuWhDeSpnRWOg70jq4eQ6E/L+nyhnMj9T71SskGzbas9HO269H1U9rrhjb19lMv35ruhPNgWy2NldLxvu4/1tUOG7zWMzw9Ke0rBfIwziSHo2d2ZtK4rWvvhMikZ5ejhbdCA8eX86zealPl0ixI195x5AU471E+BwnPy4kJ8+pOMR+Gzakz+iSnI350wejodYUUT2NhNhHbclfQ+XARwt/Tls/ZrkfVT2k5W9pUslz91Mu3smUzKccxUKdXCl8o9OOl4tlpeaiQ00PFs5pntfOku8dBzSblvafAeMWq5rnu0zsmE0S67qu7PRfGuOLKC6R8ObZEFE5rEQBLCtsgj0dnZhWsYpghYwRmFX/r2RqLHv8lxbOFwzuh5vn3f3PE3Q4IhF3boB5bP0WbTW3q7adevjbc9VTJ7yf6uh+sVLJcjD9sl9kE+PMsbdWtZOJUCOkwQ0HYM/H0dRynr3I8Esygygur4haR4plBYnhOUvAjpuLWmeLHW3RUmPdFVY+pNLwtHt08LeC+VfwAXD1TPjxvFXBuX7fkKr1JytvUOcvMzPg3BWYGzEgOuj5SsKUwUWPCgf+oMNBhzLT3s+qIjtZosCHPbB8Mtg317LLvnPpKZth1DmALeIXxeOX9UnHVvmZFiOdo+il1Vn336qvefurlM4g3vWIvBXZaGEtb406tXMbM2aSyGG+bi5GB41EGBhLerXDiwVNRPO24MSrqPUkGFUcWMl0rGvGyOqDi7FlWSXx2bJBPLwBs0fkp7WUnxF7AVWU6E5o6I0N1ocGDJY72Ez0nJ6rrVOcGSwZ43g0tMrryryLKlgCcszmeg+6xB7OnzxWSE6UQ3S+xL/KrbUwyw66AuwMJW9olEzHsa7NgngmDSaHSp+x00/upVErvMarjEwwdJug6p696+6mXj6rsQtLvu86Ws+UrX20lNSWLo9SMDVUaOB5xMUuncV5WcwwT2C+sDfg5J430RR5RuldlKZvZGAcVmhUnv3iQC/WD4fvH93+VTv1YaW25gvDoDI6Az3KT+9Sxdc+7nqpDVxorNpwT+7XNDqP0KklGsqOuS2SgH3v5OBrD1VY6z5VmcXn5Lvt2GSbbmMoNu+bodvfCZY1dkUI7ZHVj9iSOQbfUJqt2Uv5j6Kfotktflf6ufurlo6JzSHJnf+dMeeYUsZq3VJ7ibpngC7vprjgRZudzBv3XIxmDR8miAUN5Y38fk/0VH42ZmVhyZN3zQddSp7CcJjsdDbRI5aEzMUgW36cY39yr5NqA2NRZctHB6pQ7cVYRU7MIZOMolzgNZUuE/oQlhOPD/oRE0puO1j0NL4o3GzTtm+XytDGTaRim7CoLncKuGZjeW2HHqpX2+yjPo/jU3/K47r5oJ/Hf+H5K/VXPvfuqt596+QomKEdJt34AL3Pc/NiLURXp8K9GcbVHeD1OigbMwNUPZGOBSqMxsyeI47EG4xmknyBLefpBr8v/RNepAZ6sc8mrMwNkwlH14Ne9BGZftu3XKpcffbHiuRZSPcHyo7xwxTFoQdhoTJP2HWXwtLGw6wi0DR7pT81+OCrjip3UDo6ln6LK3n3V20+9fCP4P/xjZ1/GYNrKrnShQjAQBTFY0rB4Gc6S/DddaycWlJxe7FedCQwdIT93DAxiDxUoDyWRwZ4/z1x7UtoUADgpZuiJxM/gvugTEJ2IqQv1ntRZ9RjMIvU8Z6sPfe4oD2GyrKkKr01XHbAf259sG5ZWYZP27WTMaWNh17WGu5r/gaKYKIItjv2NAu81+h8FOux0LP1UavnGJ+m8tK96+6mXjzofVB/GQVb12AfCRjZOpogd/9CXx1872KU4HA9OgfcROITHug4M0SgVgC4b6QfJgoeGnmbKXRlk4ZlGTODkxWCGrbhJUh5AgpJT0zOnoRggW1tz8K+hSZ3XCO/yGqZMAmi010LCE7uBMc6WerBlMCDxeO3rbmNh1wHEWz5gK4hTqP1ugO7/UeDgzSuFpp2Ufiz9FD337qvefurlOwhf+htjI4d4Ut/XlX6II+ptpvvNSeUwfjI5TL5Az96dr0V1uZ3lwuEkZbO41u3HSnzbYlBa7xykSLrvFKJR4CCS05iQUUs2B/lGMjkJxuqMGTkv8G1rqJb3IB5WFE8VcH5e8ujslVXjM0wp69pIuLA9yIkYVmtsEf6u+zGuc+3raWNr7UpbIFR/sFsA9aTtKizobxC7GeNV60+K/zHjgc9jJ/g8tNaeS/op9drbpt5+6uIT/kw0eadJ283HYfrY4F2nnjcnlWmnSdnZ2NXpUPmLTAMUnFMgjdm8eSZmcEujA3hkv5ZCtjoB2OZxaaVPETKhtwaa7lnxsN3GDKGqi/jvKd3yo4eXJnWWbNeJE/HVXhBOYdrXVTKYHZkefbxuktNSemmrEocyZ+vvIH5mw9iRbZmPunvKm2tf6lq1CwJFps8Su2J7JiLWiVjBczrPBr9UQOHPudi1tBvAxAvc2IqziaDHTgUYi1Fr7Lm0n1KRvW3q7adePvoyhKNl9QH9R4H2a3ZJkXv8URngRR//IHRBKV2h5nG9BVPJNMA1MtDoaOw4Gd69cMoMp7MF0ShZFvYzON2zh42RGHCY3RXL6uIZgMcz+Kl6TeosmTWHMiXb0g3TyQarskqO5dDpxcovDcAm2HNVHnBFxhi714rGngRzHm77Sh4N29PGFttV8pl49Ke2VCbtg47bPMauPCdtV+mPE5eazYEF2xzE57UT7B5abE/VhTa4pJ9Sr71t6u2nXr7Ur6Tz7D7rMYKDh/Lp47SBpwr9uOrIO5vldpeDQimw96y6pwEupi4/jfmF7hkgObDAdg0NcRVJRqqvhJSWoAw+EKueYyTDnY5zHZTevXX2q5bfpc+x72QbW2NX5aUuYDee1YPj6jZXBcKfcN12pW9bHUq1Ntwm7VTKXIpbY8+SvBsWZ1jSvlrk5UPGby1BO6eliaZsxuq39wN7lWmOh62IfoarwlkJ0JFbREM1UEt8diTYlEjeVIxJrsqw5WQp71ScbZ2Y7J5fcomjbp/pvlW/Ps+Mmymdi6JUDxzuFJ6W1ztDMv6tr3Qk21rLZbMVAxnmc+3raWNr7Fqz9aXqbJgmBQp/zsGubOWU2iC7EayIzK4eOxUgLEatsWdR4IzI2TbdqZ9a26Mdtoj6srV2hVSvuTszV2RMRagMdovSyVXup/jXpt/uBABOKkyF0oHZmugdUaUQvLMNRiWW1IAlh4EsJ1OqCHLO2LhPhmjU0fZLt171TOl8pcodnmz7md5XeEYRDA4MBGPcRmy7PXJ6ZrCSVF3Am3ZBw7R6zbWvp43tYVfKpe4tOnm7ym5sj7L1aDsCB92DCxOILzNwPHbK2Ju3e9izWWCWOMumHRZ79FNvf8YutpOT1KBOCoxlU2NxpvZx3F501WSw4f0Ly6yDrp59RsAAGIINRmQ3olG/sAd4FCiHwYuz4n2a8bSuyoO8HxW4ppmb4n7W/eD3RorDgEkPXTnZdE9XyvQO/GKv0pTOpYw0pjkNhxmozT5L8naNE04cb3+sYM6b8sCbb+7l9Zpr32Ibk8yt7Pq2Agzyp2x/8nYFG2GNDTk1ZbbFyfDh2rx9Fu1Efg9J1lb29BTX4plr0736qas/Czf6HWMp45f9lobn4nvcluLHkHbr3TvXIayiLgLlHyVw5jwfkIq8NzFS9WZGRkfkR1ou8uosPgZrGg1lMCiC0UvF551cUUNSOu9YmAQsfrmnvJS56HDBsDbH9SS9aY+DgVRxNPD0odaWNl3eybYsvrBrC8gd0oT57H5KNTw23duekr+6P+8A6bWLtBXP0or8pIxsuRyl41motEtnNThm2fzimJnUC13tJFi1WPEwW2SFRhlrCEdHODdiFd3PXDvseV81teIBp7Dr6bWWSZt2bWOXfrphfz45y6x1PCzZ/6fA8vxoSA3CZq0MUqwMWN72v/+YUGSuzjgSr2Nmv730Yn+iSsNk6eMtb5jxyJ+kNz94ZQWb7Cl1WMnm7y9aGoZdW+hcQ5rsuKafUuM5Nt2jn27Sn68B+t2LXLXVRu3UOHjPwhbS5Ix+d20+UAFencWHY0s/uvRUTfzs7bIP75mhe0QGzwwEwq4zwDoSVo9NxbNLP5Xc6M+VdnK7Ej8nmhnl8zkZToDXq3M/i1Ij5EUhz0VSGgciNv03DsWCIrKFQNi1hc5xpnlsunk/jf7cbiyrHY8A5l0C7zDY3jgLmqEzWz0c64RwzsWVjOSxpcCL7e9gDLoeBMKu14P7nqU6bbppP43+PG3R1Y6HIgQ022x/6sqS9SzIqTN7zJyoYjXDdmTthT/vyGZ9P+0sQL4GJcOu1wD6zkU6bLp1P43+PGHT1e94JuRHciAQCAQCgUAgMEBgkxXPQGI8BAKBQCAQCAQCDQRuffPNN8t/QdoQHEmBQCAQCAQCgUAJgdhqK6EScYFAIBAIBAK7IRBbbbtBG4IDgUAgEAgESgiE4ymhEnGBQCAQCAQCuyFwsZvkEBwInAkCOq57V6ryU4JLBe4fKvDbtuYHYcUTFAicJQLheM7S7KH0xgjw42n+kVb6oriuX+mZbxh+tHE5IS4QOAkEwvGchBlDiT0QkAO5I7l8i/CBwqWea/8+Y/y/pfgnh6x+ggKBQKCAQDieAigRFQiAgBwNX5rgo638T5XXxJVI6eMttfT/Y0q8ERcIBAKHQzieaAWBwDQCfERyvKq5kksOiC22Jwqr/pHfFcEREQicGAJxqu3EDBrqbIuAnMnjTuLk/zgSLx95faTAKgknFBQIBAIFBMLxFECJqEAgQ+BT3b+VI6l94DVj7bfn+Egk/5Cu+m8wBpniIRA4MwRiq+3MDB7qzkaAFU+/2pEz4f0Nx6U5eICDgf5WeKQ0e9djTopDCRYHX1AgEAgIgVjxRDMIBCoIyJHgXFi1/AJL53S45TP6/FvjLxSHk+EEW/6/lnBWxP+kEBQIBAIjBGLFMwIkHgOBDAEcCPSrHEy61/WVAs4IR2MrIf6X0nPFv9EVYkV0X8+28kmR8ScQCATeIxAfCY2WEAhUEJDjsJUNW2qv9RzbZhWsIjoQmINAbLXNQSt4zw0BVjmsbO4rPJPjicMC59YCQt9dEAjHswusIfTYEZCT4f3OXQW+ufZMV75g8Hs4H6EQFAisRCAcz0oAI/vJIsDhAcje49iXC3BGBzmg+J0OQAQFAgsQCMezALTIchYI8PudP+RgxgcE7PQa32MLCgQCgQUIhONZAFpkOQsE2GrrP5PTOSAOGfCuh9VOn3YWaISSgcCGCPwfVkIoNvgOQSIAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\frac{t \\left(k^{2} K_{1}\\left(\\frac{k}{t}\\right) + 3 k t K_{0}\\left(\\frac{k}{t}\\right) + 6 t^{2} K_{1}\\left(\\frac{k}{t}\\right)\\right) e^{\\frac{k \\left(x + 1\\right)}{t}}}{k^{3}}$"
      ],
      "text/plain": [
       "                                                                 k⋅(x + 1)\n",
       "                                                                 ─────────\n",
       "  ⎛ 2        ⎛   k⎞                ⎛   k⎞      2        ⎛   k⎞⎞      t    \n",
       "t⋅⎜k ⋅besselk⎜1, ─⎟ + 3⋅k⋅t⋅besselk⎜0, ─⎟ + 6⋅t ⋅besselk⎜1, ─⎟⎟⋅ℯ         \n",
       "  ⎝          ⎝   t⎠                ⎝   t⎠               ⎝   t⎠⎠           \n",
       "──────────────────────────────────────────────────────────────────────────\n",
       "                                     3                                    \n",
       "                                    k                                     "
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "edk=sympy.simplify(-Pk+t*m*sk+(x+1)*m*nk)\n",
    "edk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Replace K<sub>0</sub> with a combination of K<sub>1</sub> and K<sub>2</sub> and \n",
    "define\n",
    "$$\n",
    "B_1 \\equiv K_1(k/t) \\quad ; \\quad B_2 \\equiv K_2(k/t)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "edk2=edk.subs(sympy.besselk(0,k/t),sympy.besselk(2,k/t)-\n",
    "             2*t/k*sympy.besselk(1,k/t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For convenience, replace with B1 and B2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMUAAAAzCAYAAAAn88z0AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAKyUlEQVR4Ae2c63EUORDH19QGAFwGvgzARIDJgEcEZzKA8if45oIMgAgoyAAuAh4ZHBlgnAH3/8nqsWZWmhnNy7uLukqrV6u71a3Wa2bn4Pfv36sCRQNTa+Dly5dvRPNQ8YMUbdUdq+5c8fcUzhzl4vdMdN8qvojRvxErLGVFAxNo4JNo/EjR0YC8ozqcZrBDqO2rFH3KVQ/9Two3QzzlXyufbFucItRWSU+pAVaIby0EX2lwvm2p71NVG+xhA9HG6ZzjheVB+oNwToJ8lVxXqZIoGphWA0cix8A/VPwQ0kozQxNTtrFCqJxBTLu/FX4qsL25q/KnirNAbaD/XXF0RVD5Z4UPwtlwzLJSZKm6IGdogAEOEDPwwrMFZf8pNOFIAxXc9woPfDq67282HJiPrjRlpRiozdIsrQENZnMIZul/lGdgh05xW/nzJgXvBBRzAOdMslLZc2JAaeiFA/lYZRzoDdgSfbZMjxg5NqA4xYZKSsEEGmBQs1ViwH7TQL2r+LZiO3h/VR6cGqj+mQLtcCDnDMqfKLgtjuLKQWio/BuF7K0VbT1sOCblk2yfJBh7xAI7ooEF7MWZgAGLE3xUeKxgq8dK5ez37yk04Z7qOPyySjzxaRwoG9T2jgJXr6wGp0q7c40Rol5p5NgEnlOMCS9evHim8HAMjdJ2nA1y9Sd7HWO33HZT4ov/CXKMoan2b4a2b2t7ANGh4L0P764tayE91eH5jxRsuWTmCIFVhmXsuXDjnhtiB2nhs7/kBoEbCx4CMUNdG4g/fbHlHNnIcwMT3eeqfDbddClBvJlFLxRv3L50tZ2q3svwUbFtq7JIqx26TY69FDG1YdX4oTg63mpOISQMyd0yJ/9WQT3uv4rZL3aC8H4J6avi8MDl2qmMgY2gXL9FBXWIiR+1QWZo24BMYM5XLN7oDiNVMihNn+jbI6Wbk0EljOrm1E3SpuKL3u4rvqiEWTgh3jeX5t/Fs3nQZja3mbtLPdwEhCf/JL6EgCbGcTcKEUToMIAYUNWgiuClitgfnqUqc8olq9OB4twZlFmfQyFPUM0BbIU4VZ2V1cQR7ty6abMpeseOQ3Re64dl1J/srYfaWPPF4iZP5Q+M+doSPmYWZ1npM3M8Fl5fZWIYwAbJZe7ql4EBxO6uL2sSv5Khi3aiZbIY5yXkAisceqt0hx4Vuuh0yT9YN55x0qaSjfd/fimwda3k7hK4rV50qsHVhrfNdeuGcBgoOqOFeOq425OFZR1pt2VSu9TWCOfCKJ28I3ygzeCbxKgR+r2KxB+HvxUiez1R1LaizqkbeHfZlG0yt0O5KyO0FwPpkskB5+07EQ+WbS0mKA1GNxVgzFUW++Avit1jeaWbgCFTs34Tlzw8NvBFH57vFIjZ27aeY4QTgxpt0cBhue6D5mQzYIxxW5nkQC63NVG6bcDV5DeaajNYN553X5tiF+zZJqMTy8vEVpBXMAAuNpbSMXrK3kkgZC7gFCiF90AYTDwhdDNXB6Ej1bfNflVz0WO/j4FXSnPjYcDAxQm5z+Z2KhvUDrrQd7IoTx8A8hwiCZ3GFs5kIBmQBwPSP1bGrwpR8LiT60Z0c2zKQEPeVhBNcNArlwZuxffyn6os+waolVmjUnwYNzi5G6fKD9lRNKims+ugCmdIbW8CNJfEkOfNwkTeFP5UnamtBMoziN8pXinUBq/y8GBZp13qhstooyyXVswVHwMTXhurk8pmBfFGhzZo6B9PdFO3TyZ/b92IFnphIAJMKoC9SnGZu/rtY1PsaHSuWgYp8aSeSxLkdH3z1cifujzxKOMj8XytAO/Zt05Iuw5EpoN9PfC2cC+Ctm1JDMOev+YQNFAZA5gZldWCK1UbTDbbgsYgSIGjrUpkp71zAk8n5UiOlnCY9WjXBPq2Un3MALx12XtVEy79Q0+8k3PLp0N+2bpR4+a1L/1gRYw9o+ljU+zSpmNVX67E4KkPttr/pTw3bbNPPOKBfOhxEVjDxTO12WBqxl2GsX0ieOYUxAxA2w6lZKINRsUBeMaBc4UzWardSnixQU85PPlzSuo8FaUpfByZ9k3+bJ+QM6aHWBlkDDZ0o4oT8QhffOPcQhlnwYq30gykPjbtM8E5OUUzSyfWiQli+LttqGSgr7VdxQT0ayRu+BxMV2JWeb3SbbPHudDb6h1Z0XADRZm2JZbZEqgMeplt//XyYfQzpRngXA6wVTGe7QSmr3VnGC9XJ/VAzlzd0Fc3QDqY9LUpdsSeXfClC2HGejc2pDOeBVVjdC5+a0+YgVkNSjFntmQGrso8nkXUMSC7wBlGSNGO+IEBL1aFKE4LA84bgLWzgYJc0LM3Lh3SAj8XyCK+xCEc+YzJaXWDdCP6zVkSJ4m9stDXpqwU2LMNqGe7tAGS56HCx42KCQtEH/7RlX1CNhUpc4pKMRKAmYPtQ1tHcRbOAl3wRAip8wSDgtkdWvcVcsEZXXI2B6EZOGrEXCYZ+Bs3MJINh0efHBKbco7WjWiyKsKD7WMT+tqUtqnJz2huPPUWb/pFOWGvwJwCg7IfZ3laKe7aO74XGgM6CmpPHUrDaNBr4jKbnytwa9LmfDRPAfTPrFJ0cD76wQBkL17VGc6csXjydPhYgYOvAf3kPbJqlVB6Et2IDrQZkJylmg4H/742ZXJqnYVFn75xyIafnXPIt7ZDiF2E2guBOR2QQlBOdWed07YvrngwC3LbErtZ6UsmC8/zzD5oZzEZiSwZcYjq6a7PrxTbKtmLg2/HDdJi+u0l2DUj2UoxRAxmDWaKfZstmHVjM+8QHU3exg9kViOcwq3ESmMDVoZcoA12LBBoYPBKAQ0ZxR7oZM1QAf9o0hseQ7O0Y3i2cz9V3rWtE9p+g3TAa+ZsHWug8qwX8byO2TLb7V+N3p+cGbNSoDceZLFHnlSxMhRONmTmU7P9Bumm9tLhiN6y2uzbKj9CHVdNR60UkJGR2N9yLffHz+JXat3ulGzFU2metk+6wm93r/tLN9op+rMqmEUDu6GBG7shZpGyaGA5DRSnWE7XhdOOaOBAn/rI/k/tjvStiFk0MEgD5UwxSG2l0T5roGyf9tm6pW+DNFCcYpDaSqN91sDYh3f7rJtr65ueH/DE+tQLcOjj1F9Or03OfWVcnGI7LZvzl9Pt7MEOS1W2TwsZj9lfgTdS+fiYvX6d4s5fLnnvy4CX9nhz114AtPISz6CB4hQzKDVGUgOa/3vwjhivVlT/r4jhqox3kuyfhAmUUjyXBsr2aS7Npuky25+lq937ZH3/ctpGptQN1EBZKQYqbkizYEvUtVJU5P2WKfWX0wqvJKbTQHGK6XTZhxLbJ7ZRvf7EJLxD4XOeSP3ltA/PgpOpgbJ9ylTYSHQOz9UqoUHPCsAHILiCrX2T1TsEZe6/Kj6/Ulxe95ay5oTiFHNqN6CtwczA5zzhPmygPA4BkK9991Z1rBCUT/GXU5EpkKOBsn3K0dY4XFYJoPbdW+VxluaNFE4CvjkLMde0vbZdwi0wQgPlhcARystpqgHNzM8H3PibbfXd3BwaBXcZDZSVYhk9w4WZnxWBj4/xbaryIE6K2EYoTrGAVeQAbJE4J2zLd28X6PXusihOsYzt2r57u5LT8CGBAluigeIUyxiCa1U++tw8KNv16tLfvV2m1zvKpTjFMoZj+1S92uGdgwM3ZwtWiapuGXEKlzYN/A9feEJTN6XxBgAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\frac{t \\left(B_{1} k + 3 B_{2} t\\right) e^{\\frac{k \\left(x + 1\\right)}{t}}}{k^{2}}$"
      ],
      "text/plain": [
       "                   k⋅(x + 1)\n",
       "                   ─────────\n",
       "                       t    \n",
       "t⋅(B₁⋅k + 3⋅B₂⋅t)⋅ℯ         \n",
       "────────────────────────────\n",
       "              2             \n",
       "             k              "
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "edk3=sympy.simplify(edk2.subs(sympy.besselk(2,k/t),sympy.Symbol('B2')).subs(sympy.besselk(1,k/t),sympy.Symbol('B1')))\n",
    "edk3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number density with antiparticles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAA8CAYAAADFV2n8AAAACXBIWXMAAA7EAAAOxAGVKw4bAAATAklEQVR4Ae2d7bHdNBPHTzIpAEIHoYOQVJCkA14qSOggTD7BtwzpAFIBgQ6AChLoAJ4KgHSQ5//T1WpkW7bll+Nj37s74ytZXq32RVqtZB3fWx8+fDg5uAZcA66BtTTw3Xff/SBavyr9ZS2ae6YjOT8Sf7/reqT8+zFeb48h+HPXgGvANVCrgehw794Uh4teoqN9quwfyuOAB8Gd7qB6/KFrwDVQqwE5nGfCfaD0i9o61wVPMv8pWb7X9XpMJne6Yxry564B18CoBuR07gsJp/NoFLkSQTQ/13WvEv0saGr/eS1h4f4o3I/ge6iOO90h7fgz14BroFYDRHg/yuEM7mnq+TNdf+kadGbgid57pX/XMjAXT23gKOHrjzYNlb3SxR51LRDlv1ad3m0Gd7q1qnQ814BroKgBORgiOyLdl0WErFC4RIN3dbEcL0J0WE+U/lZEmFAoGoMOU8/hGwePk+xzlD8Ljyh+FITHpPNGV+82w0WcrhhjVnk8KsEBEIZkic9GZ3XEFO5g5ziAKnbP4pa2umH2xCEREeJwBkE4bBcQWQ45VOitNR5orxfEx5+6XgmhN6KOvD5W2ueU2/Thv3drZHOnK8aZFT+LgrSZPdT9mCx6PjqrZwJ/L/xfs3vPrqiBC9jqRthTeiVSxLGNRrnRnARbweGq7j1dz3W1o0gcXMcpR1zwf9BldUM+0j5nAj9VgaJ4w4GD/3WJoU2drphhpnihtMhMicG9ltXIIhw649isHkQULoZiGTO417VXfeyZr0vY6gbZ84Vs/5vkHY1yYx95opQzvIwNHDaBCStf7k8x7dBSOePiF6VEpWxP4GzJP9AFnXPDWzUA77XwsxDZtujApk5XrcPIWsuGjjAbF9TIUjOrJ7bVieiAXyutXcakup4Z1MBFbHVD7MnKFf3WgkWL96UfnCgOlh8V2PIe5/tvgVj+Uo3xYavCp9Ap4K9dBJ9hYqgkTKRLwIV+GnDH7vQQgt/Ee2YPBP9G5b0b3hG3Kon0OcPXmS1UhhIxXGhX95/WEBUes99DXSYYgmKcxjlB3UMbHBT3TtciuUQPXRVlUXkOpVn9f6rPLG2dLMcnz6TEJnxDBh5sCVHGVfqDaF3ETugrynFJW61mzygLb9jpv0zQk0H1iArpX/ThN7qfveqMtOCBcTcKwqdNxjrj+1Pd/62LPdXcx/CccdoA4eTy4rgD33ld5dE1beSA7c1BWzkrypyelQ+lpYmgF1/0kY0xjg9oTArB6eohjOIIkkNUnn0WfmHReYuoMhSD8XnW5zz0uAEM4KKgooGSoQVNnGIVCJ/lxUkpv2VGyMR/i8BP8Z5ZsWPQFm7Nba8srcp0jre6wqzOM7Wfz+oUtQEdsR/ILLkGr236o/dqd9X+IHqXshOyXtpWa9rTnNZno0bsQZAtcHC8U/lLKcv0JUCQdBKtWh/AeGArgtUceYIhnC/7s0aDFP9SBOExaaSoN6+rfGcCURlbGX1+odhGTyG6Nx57UDrFTEZBR/mT2/EGB9tgWIzSWRn0KKYNKAwmpnh/9jfaM06bLgodw2nUEZ/UARqzyVVR6BBhX0V4X+hay4mNyqK2bIAwqzOhBD6V5rO6sZnSyCM4X6bCGRnReayruKdUQW71/iBeLmEnRL2ordayJ4KIFoP4Y6WNscqzGTDVgZSawJkN9udWJfBtfOe+I8kj2aDH2EmgMpzyf7HgK6U57/ipLQCemKimAPjW71O92zGHE+VoU3uGwchEXA0lqAzlEVlWOTHhWcO9Uaxw4AGoWqpcoYa/Vs+MGQpFD75xHsysRYec0ajOitaoLJEYfIVZXSkTV5i8VL+ty4jeSNDT0u0FbNm2Z6ORgZtz9IdN7YRsO7LVGvYM5pJMVWNuwLZrPsKmuQMco01/DGNRcuBc2VpgUrSVqNXHt6RxojxtsBXCNhXjnHHFKQbuFzld0cCh4yde6Ar5SFe3DcDnTfUjQTeiZz4jELwTyeLoWAL3GbQ9eFH2FAaCcQbowwZCsWzo4wGcEjDznVQvOWvlEfIrpYsMUmpMZTWyUHVoVh/jixlyUaQLAwvgHP1hazsh/l5sdWl7LugKg1XxC9VOV+ORMZFA932Bha200jgRboqGRSCN9URsZkZ04T+1UyIjnOD/Im4Jpa/MdMMEklYEd8AWsT7hg4fWc2YkOjCCwwBEcNJEb2+Vhj075fvgoR6MOVPoJ2WKJi++qEd7vDjoqw+PSSDhMXOyZ5UbSUWrQY0sNAbfaVYXP32zeokx5CFS5+qTu1RvlTK1eY7+sLWd0MVebFVlT+kdHTEOsDn9ByB4oD8z5hhvpER6wUaxzutY/kYpDouxAyD/4PiMdKvxA9Xmn6lL7mbtwp14CvLBmy5zXAXMqqKl9WnEJoGqBjMka/tuVnYKTjcvsLyEpQNg4DALoAjlUQYGeqy0MWupbAjoQP/2IYgWz2mPt4+n2AZZ7nm5xtV5CSc8OigQnLXuWW6gIBOWZ2vDoCzWmHhp6Ef3fY7MquSp6Qr9M2AvDuJ/dn9Q3UvYCZ3txVaj9pSO4JUXqO1+E4IHldOncb6MhQS6p3/kL8YIhuzFJTbjZThHs0pjgjYZy7X4ebv0TcBku7pb6a944mUb+nip6/1cstCZW5d6qk8Qx2mHkv4GSasOK3dw0HOC2ynXzTCr2mHk/CmdYqojuKs6Q4pLg1JMhrxSokSYRdgUASufg3XQf4SPgeg8OGdmSJu9c/xGPtahXulXMQ3c7GZMlgx1dtZ0RVt7gSX9YamdiPif6Wo4nArF7MVWNfZ8IHk43kSfzyEEIllBn5NjnOBA07af8jZObXxlZEKW9lIwU4Gf12/zmT9bJS9+Bpf9qzQyToS95D7/M177CuOTHPFOfmN5NYKh2cwuRWcYMBnW6oykGKivs1CVQUnHhPY7E1IpnWboeIx1JmaUMFsLn0iX2emFrl4+hQ8OztpmeTpszfGSMVlOolP17ziEd0s8lGBIVw180cBWpof8WXDYel6a6dnqKNk2r5/ysY0l/cH4m2On+5l8Uwf6Xmw1ak/pmFUkeP+RV8oLo1LQo+JeeNfzpE9vU/Fz8hYQvM8L+/KSqWpMlOqrbql4s7K+9lXeN37bvDX0f6f9VIRwRneVWnSSUFRGZZYVjZMCCaE/g2HMSCUsBiUzNQ6WpRJnhnG4Y8CApLPmszXOgU6LE2Wp1UcH55y+/UkdXeHniUqHlhJjspxUv9YYffKZrmoGa8mpwgORPhG/TUZ9bQ2Wq/4a/WG2ndQ+9mOSGF25FATZi61q7Un/J1hAVvonqzA+l1i0sZ5fEqxv4hNGQTIsHROjbRwFobG9EDs2h5VTFKQ8AxdHC4SIRfcp3Fa+SulX1bt/Y33os3dD52IZyz4UA7UX9DzwIoTSBECHBejAHVBd2oPvtnNlkA622yF2ngLTKfxcDKQnBv+i/rDEThcTfFrDNbYaxYl98qSUl8boHCfFeHim/B765DStOHaugcbLxuR0o2EfKm3vozDwbFYj+iXyCCBcnplDjqXFBOdmHa+NYEejzJHbkifQVRu8HCuBReJWL+GoDmW0+bnypXZLZdRHTotKuC/BkCwl/FAmPphIanQFvvFgeu+le64H4pWBvkZ/WGKnpeJNttVEO8Ffja1qcNA3q4oE4oUVHFtkD1LhfjLvIysm2y44m2G/Lfg2XYW2wvaCGMUZEGGyxGaPMAeW6bZERcF05JPKcFxEwXSKMXgrhGLUqfLgyEWnwZjKQztKG5vQujfA4cNHmgTsQUyRg4iXdtsTSUTpJMiHXEMwJEuxnniEJlsdJlMRLyvEHux/tnWSoZwvq3bX7A/nsFOt8JNsNcNO8FFjqxocaL0QD+3/vkDfaQcWpT5K350ySU/Fh78crK0SLzneZvmZ9jsbf3EcQd90FdoKTlc5luh0jMZMGzCyyFb3OC/2WwOeUnPGEbU3wTHyBpqr7Ugw2kuryXNdtMOREcLy/Bm4r3WRwu9JOEwWjbOIKsPZmiycTLive+iZ02vzoMcBoGs4saiTDMnSQY4FbIVQrxbY22sPtNq6a+At6g/S81p2WirLVFtNtRP81diqBoc+SWTLdoKSAAQc7OvynoI+TAARol7d0++f6mIcUM5zxg7ljBnK2Z4AcOasWvgpvNGpwr+q3v0rOrQFz1Ufp+pSOEvJHPudhZFIlHEAoKcEtz58mP1SMRGpychA/HYao1/SmSRWIz+PlCZnqDzKYD9t0PHGuqOyCM86PtEeikf2n1Se2tR9B/Sco1FMbukFYQdppEB1aXPxi7SRZjZ7HOXBAU0a5MIf7XfCmWUnhFfdUVvV4GymyBUbinL9q9S2kFakXk9K7c+2X30r0zHFF2OQSZDvZSTHe2c6qdk13qgmxtmF0xUfRANpZpRSyPcdItejBlTJIpo4b16MQJsXhaNbMcJhdiQKoY0lgJGToZcQOnjdUVvNsRM6qbFVDc6B9cv7F3snczEx5tpvA4Yfqg1WKY1xuKXTZY/1d121+6tn1YkUwf90InJiK+IfXURQLNdqYKosONHayYZOjPNvGKqGqRxH9Wvby6vtLi85LIph4iJyD/ZS+qqS2Sm2mmInmq+xVQ1OpSi7QyPKZzuktG14CWan2u/cPMJPZ1W7mdOVYThryYs6ThSMRnzn1gb0xcesCUD1qmURLs7C9r9qxIKniy7XapjcCkf6C6uFue3V2mqGnWCpxlY1OHPFu3Q9VhFMavTxi47pmfY7t/7Qi+2rp7Zup9w2GSLJF9s0dfZWamVJs686BhEB90XQM17+sZc7uKdcrOyFQxqosVW1nWioxlY1OENM7/2Z5GM1xopqD0HCJPudW7fSDQ4X6Kw4N3W60UjsbbJEPDRMkIVti7dRWCacokMVPZbRvJyrXTZHkp6MaaDSVlV2oq0aW9XgjPF9kOe8KLr4vq54qLbfRnplImJF3BnvmzpdhBUTLEP4YLrNBBQfEiplYfnFiQiiWE4u9O3VsgxNvwQ8pEJ2zHSFrWrthJQ1tqrB2bHGqllji4EV3KXH8xT7VQu3AJHxzsv6Dmx2ZKzTshe4BlwD10IDcrisXB8r5TzyjQfp4bmUwNnoj0vK2DzSLTHhZa4B18BxNSDnQlTPyZLe9xXHlW4W52wj8h6hCO50i2rxQteAa2CiBnC8/Fr0RkOMcjmb23uaw53uje4iLrxrYB0NyMmEX08qtZ/fr0P4QFQk+z2xS5Q7+G7Gne6BjOqsugZ2roFH4o9fYOJ8bhRI5o8kMN8seap858RCrgx3urk2PO8acA3M1oCcDSdzOCrF/xTDCd0k4Ogcv3Dt3VYwZbjTNU146hpwDSzWQIzyWF7fmGg3TjA43KoPVPmRscXdzAm4BlwDroF6DXikW68rx3QNuAZcA4s1cOvbb7/d5oO6i1l1Aq4B14Br4Pga8O2F49vQJXANuAYOpAHfXjiQsZxV14Br4PgacKd7fBu6BK4B18CBNHDnQLw6qzdQA/E4Dr/yAewYEgfQ+77WdoXpf10DO9WAO92dGsbZShrg/GP6+r7yfMKPfxPD91MdXAOH04BvLxzOZMdnmOhV16+6/tP114hEz4STf6uVzwj6F61GlOaP96sBd7r7tc215UxOlP8Zx89F+Y1659+ZtAQnyn3XKvNb18BhNeDbC4c13bVgnO+vFr+ub9LJObd/WokT5tN5nf+yanU8dQ3sWQMe6e7ZOteYNzlN2zIYi3STFlQHJ/25Lv8PBUkrnjmaBtzpHs1i14dfthfYZqg6hSC8e8JnP/ez2jrXR1UuyXXSgG8vXCdrHksWIt0U5cqREsE+1MUnAfkma3LG0eFShqM+xXvS8N1SpdAiCqb+U113dbEN8Y8u6vDS7pUucL7SBZDnZETiIZT6H9fAmTXgTvfMCnbyXQ3I0eFYcXoc/zrpHocL2HEwjoTZfyIgwqUcp0sdAIfKv4cxeKJnPOcj0vzLmLfc81ApdDglgRNnL9jKaZNvoBb/eaDKHVwDZ9GAby+cRa1OdEQDRKbAb3KCIa+Ujz/jjNsnGnDA4JDa9Sw60VOsj7MFcNAcR3sV7q7+3I15yvOolnLac3ANbKoBj3Q3Vbc3FjUQ9nOVx5m+M2eolBMJjZdkKhuLRKnP3jAOFKdLFJyDRcftUxC04ycgck15fhMNeKS7iZq9kZYGcLZEtDi+r+UwzTG20MZvcbgRC5on3efRLEU4+D8zPMqAL3W1ccMD/+MaOKcG3OmeU7tOu6MBOT+LSF8qT1TKvuofys92vLERXpCVIlecccO5qi3K4MP2lMk7uAY20YA73U3U7I1kGiDCBMwR2q/N2Bo4ySE+J50BJedqDt72fI0szp7olxdrtGs82XNPXQNn04A73bOp1gn3aKBvuR+Of6nOJz31eouj48TB/tRCCs5Uz83B22NwrYztjfZ+r+F56hpYXQPudFdXqRMc0QAO76XhyOG9V55jXDg/otz0zHAqUqLV0k+DaavkUPmRBacZaC9sMSh1cA1sooH/A5D8UeVg4lKcAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\frac{2 t \\left(k K_{0}\\left(\\frac{k}{t}\\right) + 2 t K_{1}\\left(\\frac{k}{t}\\right)\\right) \\sinh{\\left(\\frac{k \\left(x + 1\\right)}{t} \\right)}}{k^{2} m}$"
      ],
      "text/plain": [
       "    ⎛         ⎛   k⎞              ⎛   k⎞⎞     ⎛k⋅(x + 1)⎞\n",
       "2⋅t⋅⎜k⋅besselk⎜0, ─⎟ + 2⋅t⋅besselk⎜1, ─⎟⎟⋅sinh⎜─────────⎟\n",
       "    ⎝         ⎝   t⎠              ⎝   t⎠⎠     ⎝    t    ⎠\n",
       "─────────────────────────────────────────────────────────\n",
       "                            2                            \n",
       "                           k ⋅m                          "
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nka=sympy.simplify(sympy.diff(Pka,x))/m\n",
    "nka"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute nka/Pka:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI8AAAA8CAYAAABfPelIAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJ3ElEQVR4Ae2d23UUORCGxz4OwJgMTAZcIgAyWCACIAM4PMEbBzIAIuCSAWwEBjJgN4JlnQH7fxqVLHWre9T0jGd2pDpH6NKlklT1q3SZbnPw69evRaO+Bl68ePFGpZ8Vf+o/3b8SjfNYo/pT4bbS5yUjPCxhqo3HA+ekFuBgXw+Yh0p+UxograQGno6KpLhHKrqp+F7n0d5nNebvGuQrhXclg23gibQk5V33yrsdFc9KSuYfCqezhKyhsvrwpESM+N6K75h+r+Jv4Ek1xIx7K8WNrvl6/kjhh8KoQeCTvHPFf6XNbCandjA6ffvWbUFlrxXYx5UQXved+EeXrwYer0opipmG53npiwYj8TI7TxRw81nyir+r+EuWYWKh5IwaXs/pO2DF4ENG/yg+lqVREg+T54PC6PLVwHOhRpTK7ERxoyQeliFm+RgwkDdq8NFG+g9pc5DUl+8Kr8Uw6OV8f+8oHgJXLJ/+jy65RzF3rWkpk1mLcVZ6Ha+jO4odcFSXenitq0o/9c+JMNLjKB+SKrfl7poKnZEUu7SeDRo/CJiXoN/0f/QKgn4owMsY4nEpu6RDS1QeP9P4v0hZK72O19NdxdwBARyAxzLGXoP8wsdZWXoGcD4pxkuw9L3x6ZtKI2vTdKYG6H8JfRQTS2GWGniWasFzoKhSYuZC12V4gABQuFwzrwGIfsKQoXgDzfLx2fM8RFaGf91F9NWBvEAwnoflGf30KCxbYmAgKJAZ8FN53OhOkvrGDGUzhxI+KJ9dHko672XB6pahVXXET5voCv1cUx73zn4j3jzzHCP1SHx4KSNA6Poe11eavVLXwNw9GdCsPhvgWJ6Vj8VDoO7VkWzGxoTAU/WAHYOHwXI64Jj3VWESqd4rhezaOElQAbPawVA3FP9QjOufQ0yWhWSZ11glC4OzxD1WIM2EA0SnkQxkAaBBEi8TIHihuL7SvcmgMpbJ0uVmsF09AJSlY0UOk8rpiExMh3HGpxlUF+EZtosiBq7cqLIuuNeamqKEoYYxSOw1hvisHH7TTzyLg8GlD+R1PccCPSn86wU9UBz3/1Imnu8Xk66U4AUTPUrAo4Exk6AiF75kdf9m18To+S4nGXNsxFV9ZZI4Fy59ARKWLDaV7zsVcfkJgJSnHZZZNs0AEA/2xOdngUcyACYnNzb/Lu3lKptQdglKONKM041k9QB0lPK5tQ1XyhJWRF4onZ669hbJvwQmwFAMHo03WTqUH/oNDJ3gjRJQiD94KD2bOklVJU+SyxiStrqc4nGrg+ftPh7Km26YCImH7oKHWRgGpEbwKLcUaPSp8gmo/PMHegZxr2EnlvdKhw2W0nbzCd9V/hG9VLmTpxhUv1Ogg9xsonjzZrR/Jh6OtlnSM+oV82eETHHjmer9IvUJr3KPvimYAfqM5SXrkGGALm/1YmKddCsd2Ps8GiAAYT1mI8jvO2YM0MYmGvBkvYvKUb7bRCpOyMsBGGxwnQIUc5rg9HAjZlYeObQXwKcygEX77mSjOJCe2d4jnDrG+ENFnxAvoKNNjBzA3uWbk5dcDBYmyhxZc+qqH0xgu/ibJEp1eekL+ycT+DCSgteBAIFLK0ahgAqjB4+k9FzCQ3FHgvFioh08WDCk0uYqrX8xP2lAGEBdwB/XZ2wbJfVndCnZaOOpcPZac2xoK0aQehRSy7M8ywhG+moNeWMkHiKqszKp+gDBgUFpjHVTIdk3dIR87eQtO2Toqfwmj9hcMeNeSer/b792qbor5W+aIdcHlR0UttvTfwweQMPMBygsMVyb26wvlJ9nkxwaxn1jrPcKZwq7QD99J3qKyXVugqJz1feuzC1b3rgsIazNnAZYVngdkf3GZFI9gOJIafZOfyuwL7G9hdv7LDnav/8TDfQOFQ486vx9PwBbE20pAFALGd1+BfZsvWjM7QPE7nprywWyc3ufXgMbKrB+h/5sqJ1JYqUTJq7T/aSKm2U2XYVWDDzsQbjs6jKYh+htloKEZQLQsZcpJfNox6oQK2mqEafyd/s3adnqVt5EXjZAJ0wo0/0mmimWGYHYdBXq2p6HDr+0UkCkwCmBYzvuKjwznjiGV4FlCS/DoGN+9lDP9AzvxTPAwukIQMLPsZ0jJGlARdukkQHvYwUIGbcUWPrg4ya1iJ/KOZIcuxDdpR+B2XuuZa+ZG/NvlIEN6HwZXfwb7nkuiupKCUDcIfEWwdgJcONKUfs2UdgjYii8OfddWwWS2qc/TOYrSicAMs+jZ9US+zvb821NCTIMXhkPjufh4PJpa51JG+aGn5UhAQ4shylflTk8z7GUc7wjo2cptoPLLnSJ/mS9XwPP8rc0jMSM3yp5r2P7sK32JWocvdjPQFFx8zwL746Z6Vvd83irBK+DJ1QgvzVS+zahsp6weZ6ladgQbn3foz5w6jtbdsmdJrd9XLcrnGw/GniWluI1EGa6zTRvv0uPuLbg7QH3cpni3ib1kntEP+Jrl6T56o/qpg0Zip9U+EWfe6nqSXrgXo67tStDymiex2tGSuJSlBe3trrPGDLUFsq5hH041m4DT6odAMSLa1WT9zrc7YzeNTXwRDCRstxLZYpZ66skjf1UA8frDL2bHfTSwBNUERK3leKmFyVWRRozF6Xc6fD1avaEFSukgSfWhtJSGiccjqj80Lsrt86dXm4sy5UFH2+OLlfWegOPaSKKpTxmHW67Gu/jJwrACe+DRyrJJttRPauWVliigeZ5SrTUeLIaOHj+/PlvfxGQldgKq9FAW7aqMfX6B9qWrfXrtBqJDTzVmHr9A23gWb9Oq5HYwFONqdc/0Aae9eu0GolH1Yx0QwP1N7O8RM/fc1z5e9CGurEVsc3zzFc7bx+eKvS+qJwverclNPDMtw8/oma/a5overclNPDMtw+ep+hX6PlN7ZaE6m6YtS/B2LxqypeQvGZ5osD38P8o4EX4e8evFeB5oACR5hdn9wmKl0EdXtlAHuW8ynFGXcVVUI3gAQS87MVLT4nBVQYY3N9lVMxSZGDhe23+/6nkZXDl3Xfcikv/upbE7A9VtWzJyHgJ+/qRTS6f28SeAi8EUR5/6EY5wOoSnir7KW6XcR/zR/s4qJEx8bcWzxUAAuCxP99iVVieoO4LUXyOkwMJYKxyv4OSavM8LFMQRl90vAtFzpMAMDIR8TVp7ImoawA0Txax15GsCjyRSdkID3mSLkgAGkDha04DDckeAD2geFYF1QoeDN8FyZAnYWnjT+6xgWaps2/ak/2OnrF55nk1VB14PAAACn/SNyYHCj1PQCUGeK2MP7Nn+yE20e7nCJU54CnOebO4jb1KVwceWQ/vgBfpGhoAGDBiI/MNO6cvvt12S5d/yNellPOB4H3F8anNs+x39B+3/WeS+HAaXwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\frac{k \\tanh{\\left(\\frac{k \\left(x + 1\\right)}{t} \\right)}}{m t}$"
      ],
      "text/plain": [
       "      ⎛k⋅(x + 1)⎞\n",
       "k⋅tanh⎜─────────⎟\n",
       "      ⎝    t    ⎠\n",
       "─────────────────\n",
       "       m⋅t       "
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.simplify(nka/Pka)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Entropy with antiparticles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAA8CAYAAABW6aGvAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAYOklEQVR4Ae2d65EdNROGj10bgFki+EwGXCLAZAA4AtsZQPkX/NuCDMARgMnAEIEvGUAGwGbg73206kFnRprR3M5lTqtqVhqp1ep+JbV6pJmz996/f7/z4AhsHYHvv//+J+n4SvFvW9cV/aTnA0V/6Ppc6VvyPDgCjsBxENAcdPtTAf39ChoncQTOGoFoDK4VX4QzQmdJV5yQJ7reKo1z4sERcASOgIDmH86I258K7N0hqQDJSc4XARmDp5L+U8Vfna8W0ySXzu9U8wddL6Zx8FqOgCMwBwG3P+Psjzskc0ab1z1pBGQMPpaALMifLyWoeH6p6+FS/KbwUfvf1NYT7c+ifYDctXWczhFwBOYjoDm3SfsDMtKtygaJbpT9cYdk/rhzDqeLADsDP2tS9L5DofKnuv7U1TvJoBO/W8V/ra2y2sCJQK637baU96MutoFrA7tDL1THj25qEXM6R2A+Amdrf1Ade6FrCRtUbX/cIZk/6JzDCSKgicSOAE8oN0PiiRYv/loXRxzZwORUwReKf88SjMgUj15nQuXIjfNDmyUn4qXo2P0ZDKK7FdGvuvzoZhAtJ3AE5iOgOXey9gftJN/BbNAY++MOyfyx5xxOEwEWa3YSWIx7g2g4guFpoM/ZgF/vJO5tZL+Q9opBcrzT9aMIijsxUdZHiksOS5s/8h/9uKktlN87AhtF4JTtD5Af2gZV2Z+rjQ4GV+uCEdAizQ4DE25wdyTC9EhxcEZUl3o83Xyo9LexnIjF/1lyH5LKs2Oej5QRJp3ikFZZ0aFo85l4j8zIPvj1ELLogh4dUr1068ERcASWQkDz7GD2B5nV3snbIMlYZX/uL9UJW+cjQFmozj5sRY+Bjniu8t+l6+DuSOTzhWJ+o4Q+xphwhMPZaejzGHd4KR9D8JtidjOudf0U058qDZ+1w2s1gOy14aUIOQraRIj9sgldLkGJC+qvg9gfxowwPScbNGh/3CGpsASx0w+xwFRIM5vkYdRnNqMTZsAOB4O/NrDLQPhY2OBg4Hzwg2K2w4Fj8g8ErZC+4MrRyatY/gQ+Ldo1bpEzOE2VzNkh4WgKfM46SAcM8Vbm5Fn3xQjhL8H2AMeh7A9tnZMNGrQ/ix/ZyFBgIG1LmCdFDPm3ys++MKh8nth4C9cWhbYht8WgyEN1s0G8WSRYmIIcumcrfVRQHQYX2/c8BZ9FkKwYavT+RGkWrSbonp0DFl52ANgJmBXEg/7hSxD6ZxI/1UNeXriE16+67xyNKL8qRF7QMvgHg+hpk3HC2PhI92wt8g5HOl4p38NR9zvRpPoyfoPcaV2lee+ENtLA76KY82L5vKSa8rP8vjjnJBXpxR/dcLLYVWnPs2K9pQskA3hU24h2+6rfOydVftI2pa3P0vfSn/n0THF2Hq2JT58uandR29PX1piyiNfZ2R90lOypzThpGyRZB+3Pog6JGsTQsG3dbCMrzbk6vxaZ/UJB+QDKp5n/Kn6juPMDVspjcYUHC2y6UCi7HETLIkK7LJhvypT5EtVjIXqu+JM8xenkRlmZVCxSOGDtRbARVrS87AmeLP6dhbYhrEvQDjhNxkgy0Kf07Z+Kr3XNCei+Ey/b3RjixSTGUGLASTPWcEx4mjMexOiYDaJjAWieVNK6SncWBeVxPNTMkSzTukywNxnratw5agGj2gpL0oGN+I2yEWn7qj84J0VzsjYl1WXFNGO4aCfXxMd0iv2E3cX+NmNU6SVtjzU3K5ZMZ21/UF46nIsN4kGxaH/uz+rJbmWcjz0DLKB4EmLRY5Jkg2hsUWs/NRq9fd2wx9sKK2I6q8S7rzr6WNt9dIuVCYtHukaf86sOC+JXusDolwqB0Av9ZgW1xwD7ILY7i5cqN4ZrBiMW+qIxzvCF3sbGP0l5M9akG/wYo01QHg4LTjThsa5Udnv6D4Ur/kEmnLgxAXrmw7HCJBuRCFs1J+kf1cF5sb5NWISkzeumn9sEA/dTbUqHrWSdNOc7jJQhXvaCY664yTsAPjj39EE6p6z9RWyPMVswTufwVLYHsT8IRx/qOjcb1Gt/7k9FvVCPQcgPTLWfJu3siAGaC9QjlLbZrd5Y40unDfG+azn/92vVT7fE8lTL5oJdG79lWxC3qBf6zW5LPG4XF3A6Q/p7jGFB/3B8IT1wPDiuwSFsO3VsN9o43ClNG+wwsQCw6LHL8k28n+WQiAeGhoX3ua6Qjnx1uxcwfmOPXgI24ncsp2SqjTDFa+fk0Ly3vjy0TTE90pgxuMQ8pE+Zi1xDYTV8YsOMTeZMRxblYVMXsT1DSh6h/CD2B72E4znaoF77c7Vwh+FQ8H5CZxDGdh4U2mPwAjALQi7wFAPPscYXXvBm96AkEzSdIPovlRnA6xRuJwP9vtZ1aKdrTQQZY9X9pn4OY88E0n3nyDCW2ZN942yINn26LjnTxro6Fl/kb9rJVRRNmEuRNkdSyjNsWJBL861Ud4n8qTZiJ13HzMmTsylLgDfA47Ew4l2udFyWqqyJD22yMPfZ6y3aHvQ+mP2hsVZfn4MN6rU/izokAqdkzMPTmMpLBpDB2wFT9HTuC13E6VcPuq0Oe7zFE6P2mS54MnlvC5yYsB2ZjDbKZkcr8HuiCyP/ONK8Fk3fhIxkR43QDz2LDol0oO/AEJzAjIDh450P9OUojpgdgtD/sQ79Rj6/EMpiDu4EsAKb4kvCkW81feC6/2f0U+9+9e6dZAr6IZsum1RdwrqcufVpxRykuhb/o7K2r//LOlxK2IUxkmlxyEZQpXdOtnjuzXsrU/uM4WPZFBNj8Vh6sVNnx1A1/BfHRzLAE2cIjJn7PJxiH3LzfdD2qN5O9eHFWL/lPoY925rQ2LznBfWXym/st9JFO2ZMLRYtsm/Z/qCq2QFTe0o8xQZZu1n7czVFijF14kCgg7NPfLGcQbdTOj3/ZOGiHi/AlYwY1YpB9eDLQAwTVfc2yLjnhSuu0mL8qdVTnAs/iF94ElGM3Bg5JgpODvzprFN3SJjAGJFskB7BGChmIWiC7k1vBheOCTg2Qfc4nulLqhim4IAopj94oZbPa21wNnWVoM1HI+ibuqrDeCH8cxct+1f8efGVfr/RdTuVO3ym1qWe6uMIY3Bz+PWyVh12C6EB55MIkocxUbQRiZBDczKQRn5BP6VPyaYkqiyXlI5g17xUPcR5LXzEFwcAxx07yxzesxstuXptD7Sqj17YluYzeuVhr8LYV5rjVcYOTk/75VleHA92RzFjgXm7J4/uc/MQ2k3bH+m3K+hOUVVQ/Uk2SPV67c/qDom0Y7DYj0fllLUFEWO/Z2B1z8Dmn4LtdDWOg9IMmueRGYOWwKC9vUs2f403kySkFSMLg5i2Gg+6qfFfgjayC5vq0xk4HBZoF1mfxIxrxWl5zD65CP0Mv5xwLAB8ovpAV4ptcPCSClmcVA7GTO7GoVQaI0JV+qPpUzJioL3GcFTQWz1i+mzVIHlwOFdvZ0CJJb6O+nCgjUMWD9kIkwXcS2PNaIht3lfbFCrFfv1aSeqVvhoz3lNsCs2sEZA3+8BXaMx0qMYnYlNjc2mS+fuu0LZlD9ke6BgXvIeSPtjZ3LuNjGzs7K0dKgMPHnyw8djjGjsmskDn9gck+sNcG5S1P1f9bc4r1WBg4WJANQtShiOdn/XuGYi62Clhl+SNLhvkze4E/JRPO3jS7d8ZCbyVzwSkfnBAIp+SwRFpCAxiG/Qxq4nglU4AeLHQBnrFffoGJqJBZjMMDWMlaHen8pz3ThuDvFNmA2l0sAneIVVbGF0Mx7+kFb/SRZ8Uj1s6TMqfW5fafZPhQVaJPiUP2Ckj9ENakEtLj/e5/Jo81a0hW42m1L7y71U2WoPnTvyg+0NXFX1sm6+9bK72iiO6GhthPPrmpNEQj7YpkuNj1bP52KfrZJsSdbU2Unknz3nx5OEIDMeE0fiIea3NRQ50TJ0I8tphyPbQB/TJnq2RvvANvGOf8UD1WtdeUJk9+HC8zENErR1z+7OHZPlGmHYKlTfL/lx1OC6UIcGYKNeKGfx9YWjw2rkgdGbknopvekbIbgR5bNEZDW1Sh4GPw8ARAo5NWq7s8SHDg6eqmzGcxCPncOyUz04L7ynsTcQxvEfQ1hh4sHuuC7nAGcPE78Zk5Vf5MQNGh9C3oNxR6K90qJ08TZ1LSwijW+k85LxPgkW8a23EWP6jbYpkwS6wiDHO+8JkmyLe2TkT2xw951WPxZjdy/ThqE92KxuNjypW2VzJwtxDLh5e+sKQ7YEH4e+7KPvXaBijpYBTQziEHXP7c4f15L/3J9fsqRgnGD8u1TzNK82EswEUauveBkvf4DWHJnUkmNglT9Z428S4UTvQs7XHFp61Geh6/jC44NEbxI/JDV3w2iFWHkZisC60Rw7IaJOoI4p0CP2lmCcM+pMFHCyD89ep4BnnhIA5+keRWWOJhX/QRrSEG5yT4mvze6xNaTXVvRVv5gtzYqpN6TKdl4Msn0kuHs6aS3lg8Cjm7R0dK28qPoM2N6qCPdypneY4XOmcLey1PWJhTlZ71zs2EyKjyfE3uuZzfcnhdsxQOX6ctT9XS8ulTmfAM0naZ5oYoPY7A2HwKr8ZvKk8kRf1eHJpaJRu82GyMPBSp4VdC4LVMweGSQw/fjOibxeCwQ7tXlAdBj/ODTsF8Kbt9pHTc5W19RfZyYVrSWSTOiccfQkGDU7Si90RnMRPdaV46/bo4TZKgF4nE4QXx4kcY/RhfWh5DatDt7sTDmNsRCpfdk6mBEpPsiktHqXbuTalxHdSvnDE/ph9a3gonx/L4ri1eSBsCifiI141NpdmsA2NXVA97Df91uQpTei1PaqHTaUOdqYT4KuL42PGMW02D4QQK9/GAbb6UHbM5pTbHzqhPxhWe1T39+5m3mgQsHgxANgh4HikuZTHC1RtIR4rv72YK6sZUJxdMyg/Jy8XxJPBxqBvbyuHiZFp0xaFD3P8kjza/Sy5tyQDnesf8X5AbAXEyqPsdZp3hLTpNjQxwAw9+wLOFXqmgfu2IWzTQD/UfspzCn27vvVFTpY27UHuI3YcJdq4O0i7pUYkhznZhlWJdJX82P4YG5HKwVjNzcmUZrZNSZm10nNtSovdareM/9IcmI2P+rBkc1GIOR/GuuiQgZ3xnI2psT04VNTnaK8JumfXx3iyNvAja8iUBmj4mXqzUzV2bKy9StsjbXOqhH2bfvV76Y8sZ2N/rhZGhG1SDN7eAIpt2ADaCaRgkJQfBlG8j2Qhgged23zulRZaWvWgY+Dxfsit5ceYjrixPMp1sWuBY8R2UVNmNK34F90jZzswwHlawPHYiRf8+MSMF8t4EsZR2fPWoTtEULsmb5BNbYbP4xQjX/sJB5GgY4enFMCUek9V32hwdtgdYkeK/uP9kvAUo3vaf6KLfiGfcnAnH7zJt/YwEOyksXNgfKroxSMbxIe2kLlvmzdbd8VMMG7G/ort1LJmXhDA6RihykYUBCvNyZ36nTGGbowhuydpocqmGHEhhn9jNxhrusbYlALbZbIlCzYIPQmPdA8mrxWzMC+Cj/jAv2RzaRc8eBANawBtk5kJQ7Znp7rYmP+pLvYG3WybH/7m9LDbjXPDcYyNaWTk3pwR8mvsmNsfAbVyYA4RrK/u7uLfe+/fT/7QYI/RoW802GzQhQUu3u8UL/YkKl5MABbMgy0oaovdntEvuI3FX+2AH47KKS3eY9Xo0Esfcwp5mj1aiPgyNulPJh/G8RflH2wsqb1OUPvIw+LE/x/KGoVOpRPKkMyrzsmIDwvgweZFbHP1OT+3GyXnIjY38tmc7QFf6eb2p2egCZ9e+3PVU/dki+KAxmPGCw5PREpj/PHOlww8CcA3OD1LMu7hxSJxiIUCrNBva4F3heys/2i6aVziGDM+eRLkJcij7JplAODIgyfPQ4yxTPOzs44xJ2cLPcDgUHN+QIxyscYLzshSNnertgcA3f6UhxElvfbnfn/dky3FC8XQE9vFscKiRlb82ObjyYXJeJCgtngZLXe8slj7UR/0WrWdxQQex4jx8EC6PRhXbTVqHGbbOl6tkRGMkeeouzQjZO2QxjG7+JxkTujC2XmuK6R1/01HgBUy1M7qc34BsRexueAc8d2i7QFmtz/9g63X/pzrDskH/TovWsqLVWxxH/UIYFGN7p50Drnrs7D4vex+VSlPcjisR92VkPFFBnuvRcmTCMh07n2/+JxUX4UdrZPooRMUQvgsZXOZm+c+/vp6yO1PHzoD7y2e6w5Jv8oLlmoisuvCi6sHeVpaUPQsq6gH+iz2rk22oSNlxv5iR+IUHMhmd0RysWvD/dGC2scZIZzSjs2dRCP+xj7ezJwcofpZk6rfsKGbtT10Thybbn8yI1XYDNqfq0w9z2ohICBZvEtvi7eoT/tWumxCjwGU2dGydw0GSFct5sXI17EFjgJuVm1tmDlOGl8lnL0zuqU5Odxt26C4ENtDZ7n9yQ/ZQfvjOyR54Dz3vBFg25QdCfPIj6UN29P8GimfQPKFzaLvOE1QCjmO7RRNENurOAJnhYDbn3x3Ddqfs/3sN6+v5zoCdwho8WeHhN9i4DcKLj4IB7bL+e2Xpd4FuHhMHQBHoISA2599ZGrtj++Q7OPmdxtBQBOATwv5WuKo722cEJwcGT05IXlcFEdgswi4/el0bZX9cYekg5tnbAgBnJIXG9Jnkirx6YTfHjnqV0eThPdKjsD5IuD2R303xv64Q3K+g90lH0BAEyH81oHi8DPWA+SbLJbuD6UYTyd8KuvBEXAEDoSA25/gjIyyP+6QHGhwejNHQ+BztcwvpjIxLipIZ34cjv8dw/+EOvsvay6q81zZrSDg9meE/XGHZCvD3vXIIqCFmC9b+NyMfzTIAn1JIXx+KL39qOaSet11PRkE3P6Ef4xYbX/cITmZoeuCrIWAjAK7AxxZXMwuSXS++CdxW/2J7rWGi/N1BBZFwO1PPZz+2W89Vk7pCDgCjoAj4Ag4AishcO+7775jG/sPXWO2s7+S13e2/6BrJSydrSPgCDgCjoAj4AhMRGCxHRI5KO8nyuDVHAFHwBFwBBwBR+DCEVjMIblwHF19R8ARcAQcAUfAEZiBgL/UOgM8r+oIOAKOgCPgCDgCyyDgDskyODoXR8ARcAQcAUfAEZiBwNWMul7VEXAEHIEGAb1Hxovx/CoswT6x5kfZbu+y/K8j4Ag4AmUE3CEpY+MljoAjMA4BfvfkmVVR+iel3+r6yPI8dgQcAUeghIAf2ZSQ8XxHwBEICLDzoeuVrn91/dkDy1OVP0rKf1Da/+NyAognHQFHoIyAOyRlbLzEEXAEhICcjFtd/Pw+v3j7ew8o7I686Sn3IkfAEXAEigj4kU0RGi9wBByBFgIf6/6mldfcymlp/0w9DspfyvcfUWxQ8oQj4AiUEPAdkhIynu8IOAINAnIq7Cimb4ckpcd5+VLXJ02mJxwBR8AR6EHAHZIecLzIEXAEGgQ4suHoZvCLGdE8FC3vj3xSQ9+04AlHwBG4aAT8yOaiu9+VdwSqEWCHpNkdkaPB7sdnuvjU91tzPKIzwj0OzC7eE/P+CTF82D2h7hNd17o42vlbF3V4efZHXdA81kUgzRc8Tfsh1/84Ao7AphBwh2RT3enKOALLIyBHAKcDp4DPeHe6xxkh2Ge9fNr7s/LZGSEPhwR6As7GtyF19+cLlVH+SrcvdL3mniLFvIPClzzswvDuieXT3ktdH+jy4Ag4AhtFwI9sNtqxrpYjsCAC7GoQfpeTENKKf9M9jkr65Q2OCeXEdj2NDsYu1sURIeC88Dnxj+Hu7s91TJOf7oaQT1seHAFHYMMI+A7JhjvXVXMEFkIgvD8iXjgbb8xZUMzXM81Lq7of2sGgLu+h4FzgkLB7kgbbVWl/rUMb/qVOipSnHYENIuA7JBvsVFfJEVgYARwRdkJwDJ7JoTDHYVQzOCOxAvx2uk93QcjC8XmX0JFH+FpXmzYU+B9HwBHYDgLukGynL10TR2BxBOQc2G7GjdLsaPAux1ulJzklUUBeVs3teOCo7Dkeaoc8ZLD3V0h7cAQcgQ0i4A7JBjvVVXIEFkSA3QmCOQr2S6wcuezkMHxDPDLkHA9zfOwdE2OJE8SuCS+50qbJY+UeOwKOwEYQcIdkIx3pajgCKyFQOkYJn/GqzQ/HtBudCpyPX1r1gqOhcnN8rBhay+O4qP1+idF57Ag4AmeOgDskZ96BLr4jsDICOAQ31oYcglul+RwX54DdkabMaAZidjlyPydPOzlngx9Y46sb2grHNoo9OAKOwAYR+D8yLNFEDSdb6AAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle - \\frac{2 B_{2} k \\left(x + 1\\right) \\sinh{\\left(\\frac{k \\left(x + 1\\right)}{t} \\right)} - 2 \\left(B_{1} k + 4 B_{2} t\\right) \\cosh{\\left(\\frac{k \\left(x + 1\\right)}{t} \\right)}}{k^{2} m}$"
      ],
      "text/plain": [
       " ⎛                   ⎛k⋅(x + 1)⎞                         ⎛k⋅(x + 1)⎞⎞ \n",
       "-⎜2⋅B₂⋅k⋅(x + 1)⋅sinh⎜─────────⎟ - 2⋅(B₁⋅k + 4⋅B₂⋅t)⋅cosh⎜─────────⎟⎟ \n",
       " ⎝                   ⎝    t    ⎠                         ⎝    t    ⎠⎠ \n",
       "──────────────────────────────────────────────────────────────────────\n",
       "                                  2                                   \n",
       "                                 k ⋅m                                 "
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ska=sympy.simplify(sympy.diff(Pka,t))/m\n",
    "ska2=ska.subs(sympy.besselk(3,k/t),sympy.besselk(1,k/t)+\n",
    "             4*t/k*sympy.besselk(2,k/t))\n",
    "ska3=sympy.simplify(ska2.subs(sympy.besselk(2,k/t),sympy.Symbol('B2')).subs(sympy.besselk(1,k/t),sympy.Symbol('B1')))\n",
    "ska3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Energy density with antiparticles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzIAAAA8CAYAAAC914/AAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAab0lEQVR4Ae2d6bUctRaFL14OwJgIHmTAEAEmA4YIbDKAxS/45wUZgCN4QAZABBgygAwe3Az89lddp6waW+quqbu31qpWlXQ0bZ0t6aiGfuPVq1d3W7pvvvnme5X/i/yft6zHLZZ9a9irvY/Uz7/p+FDn97fY527z/AiYR/Nj6hyvAwFz4zr60a2YDwFzYhpL4VO8TnswneWysXWHPpZvI2ZZqHu53yL2ajPGy1Mdf+gcstgZgbMQMI/Mo7MU6IoTmxvmxhWr90lNMyeOc0IYFa/TNjNkVNln0oT35X9ykkY40ckI3DL2avufAu5bHS9OBtAJjYAQMI/MIxNhGAFzw9wY1ozbDTUn8jkhrIrWaZsYMqrku1JnFpMf3q5ab9PyJbBXnh/reHubFr0uVXX44vXV+JnkflDsI+o9LuUYIzCOgHRn9jEMfdRhHo3D7pgLQOBauaF2Zc0vdJFkPcdcgK6uVUVzopwTJRzaxJCR8rAb/oMqyi0ku3URyMJeffNMx186Jgdv5FT9e/l/r9EMlYMBQt3+6JansO908M5VjuNO4AvJ+xGzHLQs00XAPDogYh51NcPXF8sN5gMdc8wvaIG5YS4EAubEAYlSTmTJr27IaJBgF5zdzOfRw/bXQaAEe8myo/RYB7f4Bp1kMAI+kv/roEBhoPKZNEIUj95gOFHumAHyk+S42zfpJHMvgR91MMDYGYFsBKQ72WOYZM2jbGQteOkI7Jkbqttq8wv9qPI8x1y6Qs9Qf3PiNYilnMiVX92QUZNYZLJzDsnt1kUgG3v1D4+4sDs1ZaSQ3+TkUNg8yhx1qsufOr6TwOjdn7q+T+SPGTpp/tR/F4/zpJXy+e4RMI/aXWQetfG45as9c2Pt+QU9MDdumQ2HtpsTbR0o5cRR+Yft/Je90uKSHXUGE9+NWRbqXu4nYP9EmVRGjNLSZ+xCv6XzL5PMMRg+T66bU4XHI2nvKLBSRPnVueJGDZEmg/NOqDf1n/waHvXQgSxtSNulSzsj0EdA+lI6hplHfRgdcoUIrMkNlbX7+YUuVj09x1yhruc2yZzoI1XKiRz5B/1iFg35Srn/qor5bsyiMA9mXor9R8qF//fBiGHxxiMyPDvM9V3tD/aj4phkfpbP3ZPHOr6vz9/XOXkt7X5XAdQ/x/0kIR5XszMCOQiYR8MomUfDuNxS6CrcuLD5hf43N26JBe22mhNtPOKqlBOT8msbMuzqUyG79REoxZ6dZNy7mjgwSjBa+CPJuJuCQfMPAgMuffmfR7x+qWWekteA/NxB1LUyuDIy5o4Mj9CBj50ROIaAeTSMkHk0jMstha7FjUuaX+h/c+OWWNBuqznRxiOuSjkxKf8wctVCjoVfPF7zvs5ZpH6p8NGXvSNtjq98YieeCh11kmeXnC8WxIK6uwCOhXRxHZU3i2sMqqqduuaRp924GiteQqeNP+p68PGt3ArX+SGeiz3lghG4vKP03B7n/ZRUF4jHYOg5yXH3Jhz9V9U/Ta/z7xVOOanjf4XC6IlwXt5P84vwKX/MwOqlUd60DeOMOzhdHevJLxmgeoBHNgclf/UcqTHhC3XwvFQPqu5SOsaes/lU50Oe5lGF7OsfYWMevYajdSZsZp9vlCdjxS54QWNVn9L5/eQ5RmWl48Cu55cam71xI9YT6CX98K0wHRzTFL7YHAM2e3K1Dp89T0Sb6vy4HMQ25MKXvDkRYHR8YTPJocqQqQHk8Z/mcRyd814D/4De+yqVwiAAgyhxsUPfKbp3idFwlysvOQYrPtH8r/yX8jFqWk5hGCPU8T0d6SK7Jde9kCwLcOpOG15247e+rttCm/5SXR7PUJ8i7FUekwOPAH6ug3NwxqB5W0f0Nz56MOoky+TW7J6l6XUeg2mTXmE8ytboYBNRfsKAEPXMSc1AU2GUI7yEDNgo32wOUgeluQWOxOD+3qm4CyfGhjn4ZB5Nd4J5NICP9G+J+WZPvKDVq3NDuF7K/AI+e+AG8zUbQs3cq/OPFca8+4mO3kaewhabYwBF+VOn0rUkSWd3qstc80TUzZwIJIb9Uk6Myj+o88doaZSbMHUqO8MMwEOPgrG4ZSDN3vmWLAvUbGNDstQhBuvuLj3ROHb1ca26H4KyfhkIx/LOymBhoZLF+FRVSrFHPnBJ+7jBWX1DX9I/LUef6cD4xH2mI21D3G2oIhf8oV4YgbkOWXRhS1fKwaqu184RtY/B6035je6d0UmpLp6SjXk0jZp5NI3PbPPNznhBq1fhBuOdjkubX8BnFm6o7U90nPpOJ+l4zxXjJRzjK453OQad5Jdch52ylhys54yB584TURVzIpAY9ks5MSr/oM4fZeLPD7GOU4eS8/4Aipw6OohbPRg6uY4yShWENLgg2+Hq9W/UiwYWOdX9WN5F+e1cuBR79KDanRFOGCw8VsYg+N9OO9GB6IMqStf0MY/D8cI/xhB3dr6or88yZJQHkxgLfgbd6rzOV5cth372dpdaEu2LSi+V15bGDH1UwsFowTE9jv65WI6oX0rGmcBlCd88mkbVPBrBRzp8jKcjKceDd8QLKrkKN9TmS5xfwGcubjA3c5zimMsZS5vxNFOHjunuyXOM6nLKWvKUtm+RxpyYRr2UE6PyD+tyMBR4qbtR8E75XeLQQSULRbIjj6oiXGQ6lPxO9YKAQ45dWupcWhfyIm8eexprMzLX4oqwFyYV7tF4Xfce66vj4i5Cy0CRfLp7PmaERvbZvvJFf1pldRNLhrbe1bLd6LHr0EsG5DFdG0s7V3gpB6NccySQWN43j6YxNo/G8bn2+WY1bmhsv7T5Ba3YnBvCjTnmzVRFFRZ3Z+LpljQ6zpecY05ZS0a99u6bE9M9VMqJUfnKkJEyjy1Uqx1qxbMjj8IxgNA5LPgwfHjs7Hf5fGY3x5XuClNmbyGs8qjDCx346Ze0dJntWnkrTwj9gQ7y5DnS+6Gc6rJZwKfxYNAYU4lMtJeX5nlpvWmLzsGWOpAPZeI+U3jvXQCFgXcMONSxBHPyxUVdDlcz/NIeHTxby92RULJzcp4jjzCuSuoR5T4uSTSnLDiO5NdwcCS+pccho/zQqdU4UpcXjzygo091oLef6cC1OHIIqgzOUR4oT9IzxuBXukY6hZOGthH+ow76PJsfdb7Z8so7deZRikb73Dxq45FetXgqHZycb2odH5wfav3dGy9o6965EfqZ9kvp+SnzC2VE2ZvNMd2GSo/Qr6o9Ok8/oNAVbeluRCrNSXNMXW7RWrIui7qyXgrXmlMSmdBDr7sCqXE/9HJc4njMWpyIuvY49HCsjlIKFgssFKodcF2zCGcxwQD8RH5lpev8qJMs+eD+OXjHf5WG8qsFvs7jz69IyCKJ/Hgxemzxh9yoUzryJf9qF0LXsajhmhfPOHrElhzlEtd8RlhhkBwjhcesMPjIl0mm9SEEhfNCHcbfdzoony+FtDDUdbrTJJHKIQvelbEon/z5wAGfRI6OPUgO/EqmGPuBbEaDlD8fBKAtz3XcjwpmRJBXhtioiNKzkK76YlRoIIJ66yAGrHfjVCf6mv4bvAtVx1d11vnWHEEHqv6r64KhwUTDpgC8YrBrjH2d3ymcuo/yQPHoNy/pw7nG6Zq7ZunL+xWvEFDcMX5QZjGflK95BMATThiZRwP4CBd0Dr3Mmm9q+YvgBc1VfS+CG6rnJvNLjdFuuCEc0EXWLaylGEtf6hh0teysc4zyLFpLSh798ror6SFhMsu6i3ySbItPlf6kNRcFKW0RJ6bkRw0ZlcNiPP7UkHLDsfguffymIkJkkOlDNBwd1lqw6xrD44X8Ox0tg0PXlPVpna53d0PhuMgbw6w6l09bIThlNXdOEE4cmGCwpAuyaNt9LRe4teqsOBajGCDk/VgHnxp+pCPSKaj5eAHn4ZBrDB6dYywRR71bbSdwwEX9BqLmCVJ9WKwuXk5GbXk3J8UzI0lL5K3W1fYXoUtjdzxDj7M5UvfTV3XTYgGCYd7FLfI+yhGlZTDDUAlHXnD0aR2AvqfxdXD1paMcHoxtgMAxjJJmQ0Pnx/hxKp8W12/V3TwKzZjX35JHtCSbS3Wz35d/KbygyrfCjXPnF7DafI5hjFQ9qjWczhmnWZcMfrVMcaG7S8wxuWtJr7vUEV23k/liC070OPSwCw7XAoidIxbszQIhkUOx04V8EjV6ykIG110sHUKHf1FyLLauQUD9MDrYTeCuzEsdQUoMkSDe1OBa5V3Lkr4yXOp8Bo0fxZEf+bcWlQoHiwoPnRPP4vB3HS1H3joI4/ExFiwszv6VT9m/6KBNrbwVhnt58Hq/U+1LhYuwVx1epYlLzpW2RHwR2aE6KOyNzMKyMFV+yP2mI0u+LpuJotLTnLpIdoqDkUUxR5SwuXNCJnU57HZ1/0uphCNwKOUpHELf7+syhsYRysZIyuUBWQ25Un6UykeZ5pF51Mw1UopcHqE/JVy6NF7QvpvhhsYr2ttzCp97jmH8j7VMWl6Ftcob2klnzB0ca9MM0nPJs+5gnOZJBr4OWY3ZicyScwztq9ZOSXmtU9WHOfZdHa21kcJJV6XVOfFedwmELZzwHyxW4bNyIikEnWi5h60rXahwdlcfy0eBW05hZIDCsPAucSxWcL0KHIIHf48peTwHiVy1QFT98CEzuwxTjjQsvFhw8YgKBlGVx0Qi2o3738Eb/A2Z+8HYQyCkw1E2O+PUld1qJkb+N2dogFL0ya4Ie5Wfq3wnV+jSEwoj+nfQ4J2jbcp/lIOd/Is5ovTPlD+TVtx1RPcI4/GslAPZHOmko4rcEX3OSYZbiwcZVZkUMY8m4SmPlN7cAo8AJptLCYqXwguqbG4kHTfHqbgxuA5QOOsF3kltLexzylSaau0hPx3nScrmDjo6NJ8MhZEmXG8dpgjmk8k5RvGsBXPWkrGm8rorELffQ+BBGiLlgiT88WFj1esc0oQyodR3uo5FEOco5KxOecZif8pgCkOrS8rJutT1pT3Pdc5gwW1Lbq1GmWPpY8e5u3OdyofMFCbNJ4tVJndmwBvjgbowAByrR1qez+dFIAbleXMtyE39f4yDVW6JnpRyBD1j4hp1yjsmmVKO3CktYwTpq90yCiE/Ds5Tp7BqXJFvHqTAXP65eVT3Ya33RfONeXH5BJhowZbc4M47a53eWDxUX8nFWmT2OUbl5a4lY03ldddQJ11HWCknevIPAodaaT+Q332xmIVV7LhgPDSGg2Sr3YHIY8K/r+MeT8ikUZWSK6AxmNLIuq6Uzd2XQZlUvnPObjEu0sWiLhZV6UvTB0n9qhzaQNvfbwKTE8V/rIN45MLIaiQUF23CcGKAYNe9cYrnfRcWf4P5N4LlJ6XYl5dwQgq1lwG1wvyE5EslCayWyn8yX+GBXhzjYOQR+hR6HOGVX+fV44jCueuXthPDpvpQRZJBNkeUF0YKH7KI+pBf95HQrzplRlFr8iDKPNUPzHLHsFPLKUonXM2jDmLCZC88ombZXEqacUm8oNrmRtJ5R04DqyNii0RTdjxKlhYQa47uXBJjeje8Slvz7NQ5JmstWc8bXnelvXXG+ZXMFz0OVYaMGseCkgU2ixIes2oOhfGSVyRkEq8sZIU90jl3axrDRtdj7p86gjQ57jMJdRdDVTqVB7l+00G5H1aBZT8VgZI2Reqw/N+KgAGfO1W0uWuE8HhO4ECdPpUMk1HqkOGLZTEosLjr4sF1xJMWvM91pdifW97R9HW7eZQpMD+aZkkB1SMMqsBqyeIG867rkMPBSH82R1QmOspE1H1MroQj8JHjH+WH/rYwVBhxv+sYczk8IG2XK4SV8qNUnjLCRbuG6hEyq/o13uZRgrowyZ3LItWSPKKMEi5FnfAvhRfU1dwAhQlX6yUSgdWE9GJRbFK37q6oXoz/jGnpOi8qsCQ3GItz15Jed0WPnOFf+nwxxaGHNS4oNxNAa4Fex/1Z+3gQASOnkpOf9Zym5DBK7pV26vbgnWSqhZzkWGDFNafhqCMDQfP544go8CHt85CnXjpoF0TmllUTFzLhK56d6//omvdZeBkvbnGBSZCSu0QsDHlchjbjqDfXYaQQzh2YZwqTVzkMKPKlDNrP+zP41A9cqBf5sOONY6Jj9755DPAQ3P5VfBb27VSLX7G4TfVq8QKPFIBO4KK/Dlfr/mZxsNYF6nsWR5QPuoRxzTti3XaTf8ODWofGOIJOo8v06Z1k4RF3aOAHjzJg4DSPmSGTOMrN4UG1Y6h84MFTHdQ9mx9KdzaflId5JNCPOPQG19WnQ+g6v3viES0u4VIgdDG8oMLmRnTbpL85N9RP3I1/ooOxORxjKX8VEWuTO53Psg5TPlNzTPZaUvl43RW9dZ5/6euuUQ698erVyR+pKoJUyhiLGnaoFnUqi10GjIJJw2nRSuwo8zWxn2q26sHAhiFG/zBZM3j+V+GbGjUqn/oweA99tUXB1+XqfsCwrozi+vpO/i7ukO0VbeGz2hg2hUHdX+bRFEgrxNX9YB4Ja3NjWuGEzyxzTJ3PSS/7T9dw3ljVk7n+prlhTkzrVCknpuQfThc1ayzvosTzwrNm7MyOIrAL7KWILJQZ3NgZeC5/bKf+aINmFvhA+bHrg3F11U5tZIJhR45+qO7q6JxFMTtkdtMImEfT+JhHt8sjc2MdbjBH7Xqe8hzTKII50UAxeFI6X4zKPxjMfplAdjMfSckfLZN9tSvETgWPy/DYSXWu68GX95eqw07zXRz7wnazgG5uZRemXUKc+mx6V2iJRo3kiS5gSOLHwSOOu54cR9qydrB5NI24eXS7PDI3VuCGxmn+e4vHcffs0AXPMYf5ddE1b6ESXPq6a3R+WfOOzI8CnZ1gFHyRnXgRvNrxL+zcWxBfHPtcENVH9H+8b5CbbGk56lQ9ZrV0QVvnL/zf3LoOF1y+eTTdeebRND7XHGtuTPeuuTGNzzXGmhPTvVrKiVH51e7IaAHFji+78Iu/IzON3e3F7gz7ZldA9WK3guvNnMqHHLg93SE61Mi/u0LAPBrvDvNoHJtbiDE3xnvZ3BjH5ppjzInx3i3lxDH51QyZukm8UO33ZMb7d8mYvWDPBxh+rxvKI4Bbv2Aen0fduh5L9r3zng8B82gYS/NoGJdbCjU3hnvb3BjG5RZCzYnhXi7lxKT82oYMt9rYhY9d8OEmOnQJBPaCPY8XviMd4BPefLFs63czqMfzJQB3nleJgHk03K3m0TAutxRqbgz3trkxjMsthJoTw71cyolJ+dU+vxxt0cKVl/H5lnn3T/hCxP5CCBj7NrDCgw9B8H88fm+kDY2vJhAwj9rgmEdtPG75ytxo97650cbjFq/MiXavl3IiR37tOzJ3qhSfeeWLYpu+G9GG9jaujH2vn3m07Wkv1AFGYAIB86gHjnnUg+Q2A8yNXr+bGz1IbivAnOj1dyknjsqvbsjUTcKYedFrngPWQMDYC2UNLtyN4b9jFvmC3hod6TI2RcA8Evzm0aY6uNfCzQ1zY6+6uVW9zIkTOJE7v2xiyKhy1XfQ5fPcm92KCBj7w/8NCXKs/E9WhN5FXREC5pF5dEXqPGtTzA1zY1aFuoLMzIlyTgizt9X1Weu0TQyZWi8/lM+/i1NZu3URuFnspW/8IesvOp7q3F8qW1fvrq0088g8ujadnqs95oa5MZcuXUs+5kQmJ0rXaZsZMqooX6vik2o/1ZW+FmXdfTtuHHs+h/itMPAjZbvX1H1X0Dwyj/atodvVztwwN7bTvn2WbE4UcaJonbaZIYOqqWPZEefxHt+VAZAV3S1irzZzNwYjpnq0cUW4XdSVImAeXWnHullnI2BunA2hM7gyBMyJ4x0qjIrXaat/fvl4MyxhBIyAETACRsAIGAEjYASMgBGYRmDTOzLTVXOsETACRsAIGAEjYASMgBEwAkZgGIE3vv7661fDUQ41AkbACBgBI2AEjIARMAJGwAjsEwE/WrbPfnGtjIARMAJGwAgYASNgBIyAEZhAwI+WTYDjKCNgBIyAETACRsAIGAEjYAT2iYANmX32i2tlBIyAETACRsAIGAEjYASMwAQCDyfiHGUEjIARMAJGYBYE6s9q8k/NuPjkPn9My3+K2RkBI2AEjIARKEbAhkwxZE5gBIyAETACJyDAfzh9Hul0/r3O/9DxToTZNwJGwAgYASNQgoAfLStBy7JGwAgYASPQQoA7LTp+0fGvjr9ake2LZ4p/kgR9q/O3FfZuEuZTI2AEjIARMALZCNiQyYbKgkbACBgBI9BFQIbIvY6PFP63jl+78ck1d2NeJtc+NQJGwAgYASNwFgJ+tOws+JzYCBgBI2AEagS4s/J8DA0ZOz904jBs/lb4n51wXxoBI2AEjIARyELAd2SyYLKQETACRsAIjCEgYyQeGZu6I9MklzxGz8c63msCfWIEjIARMAJGoBABGzKFgFncCBgBI2AEegjwaBmPmB39Aplk3pYs78e8lyPfK8kBRsAIGAEjYARqBPxomVXBCBgBI2AEzkWAOzLN3RgZKNxt+UDHIx1fhsFSGzFcY/jc1df4vF9jZwSMgBEwAkagCAEbMkVwWdgIGAEjYARSBGSEYKzwqBifU8YowYjBxeeV+cTyDwrnTgxhGDLxpTLek/lSh50RMAJGwAgYgWIE/GhZMWROYASMgBEwAgkCzfsxMlCqc/k/Kx4DJ/2SGQYN8fhxPJPs0cfRJG9nBIyAETACRqCHwBuvXr3qBTrACBgBI2AEjEAOAjJEuMvyqQ7urLzUtb9ClgOcZYyAETACRuBsBHxH5mwInYERMAJG4KYR4C4Ld174AtnnMmTisbGbBsWNNwJGwAgYgeURsCGzPMYuwQgYASNwlQjIaOHxMd59ea5z3nf5SccfNmaEgp0RMAJGwAgsjoANmcUhdgFGwAgYgatFgEfKcPHFspeHy8q4uZNB80V9bc8IGAEjYASMwOwI2JCZHVJnaASMgBG4GQT4jPKfMli6L+zH55Tfuhkk3FAjYASMgBFYHQEbMqtD7gKNgBEwAleDAI+WPY/W1AYNL/3zrgx3Y5q4kLFvBIyAETACRmAuBP4PyIHrMvLauKYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\frac{2 t \\left(B_{1} k \\cosh{\\left(\\frac{k \\left(x + 1\\right)}{t} \\right)} - 2 B_{2} k x \\sinh{\\left(\\frac{k \\left(x + 1\\right)}{t} \\right)} - 2 B_{2} k \\sinh{\\left(\\frac{k \\left(x + 1\\right)}{t} \\right)} + 3 B_{2} t \\cosh{\\left(\\frac{k \\left(x + 1\\right)}{t} \\right)}\\right)}{k^{2}}$"
      ],
      "text/plain": [
       "    ⎛         ⎛k⋅(x + 1)⎞                ⎛k⋅(x + 1)⎞              ⎛k⋅(x + 1)⎞ \n",
       "2⋅t⋅⎜B₁⋅k⋅cosh⎜─────────⎟ - 2⋅B₂⋅k⋅x⋅sinh⎜─────────⎟ - 2⋅B₂⋅k⋅sinh⎜─────────⎟ \n",
       "    ⎝         ⎝    t    ⎠                ⎝    t    ⎠              ⎝    t    ⎠ \n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                    2                         \n",
       "                                                   k                          \n",
       "\n",
       "             ⎛k⋅(x + 1)⎞⎞\n",
       "+ 3⋅B₂⋅t⋅cosh⎜─────────⎟⎟\n",
       "             ⎝    t    ⎠⎠\n",
       "─────────────────────────\n",
       "                         \n",
       "                         "
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "edka=sympy.simplify(-Pka+t*ska*m-(x+1)*m*nka)\n",
    "edka2=edka.subs(sympy.besselk(0,k/t),sympy.besselk(2,k/t)-\n",
    "             2*t/k*sympy.besselk(1,k/t))\n",
    "edka3=sympy.simplify(edka2.subs(sympy.besselk(2,k/t),sympy.Symbol('B2')).subs(sympy.besselk(1,k/t),sympy.Symbol('B1')))\n",
    "sympy.simplify(edka3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing numerical values for nondegenerate expansion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Electrons at T=0.1 MeV:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "tval=0.1/0.511"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "Tval=0.1/197.33"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "mval=0.511/197.33"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "muval=8.0e-12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "psival=(muval-mval)/Tval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "xval=tval*psival"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct an array of pressure, density, entropy, and energy density for a range of value for k:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr=[Pk,nk,sk,edk,Pka,nka,ska,edka]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr2=[arr[i].subs('t',tval).subs('x',xval).subs('m',mval) \n",
    "      for i in range(0,8)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.78966486e-04 1.62882635e-07 3.37209953e-10 9.63299320e-13\n",
      "  3.27050929e-15]\n",
      " [3.53154567e-01 6.42832607e-04 1.99624920e-06 7.60351419e-09\n",
      "  3.22684799e-11]\n",
      " [2.79478879e+00 4.14195654e-03 1.19406589e-05 4.37722408e-08\n",
      "  1.81481669e-10]\n",
      " [1.23733555e-03 1.93611729e-06 5.71390187e-09 2.12189542e-11\n",
      "  8.86981063e-14]\n",
      " [3.57932966e-04 3.25765259e-07 6.74419874e-10 1.92659852e-12\n",
      "  6.54101806e-15]\n",
      " [1.11500790e-08 4.05920517e-11 1.89081532e-13 9.60256926e-16\n",
      "  5.09403121e-18]\n",
      " [5.58957751e+00 8.28391284e-03 2.38813168e-05 8.75444763e-08\n",
      "  3.62963310e-10]\n",
      " [2.47467105e-03 3.87223445e-06 1.14278032e-08 4.24379057e-11\n",
      "  1.77396199e-13]]\n"
     ]
    }
   ],
   "source": [
    "arr3=numpy.zeros( (8, 5) )\n",
    "for i in range(0,8):\n",
    "    for j in range(0,5):\n",
    "        arr3[i][j]=sympy.N(arr2[i].subs('k',j+1))\n",
    "print(arr3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
